MATLAB程序设计与应用-5

1、第五讲 MATLAB数值计算 符号计算 授课教师：田 鹏,MATLAB程序设计与应用,数值计算+符号计算,5.1 函数、极限和导数 5.2 空间解析几何 5.3 数列和级数 5.4 数值方法和符号积分 5.5 线形代数,参考书目： Shoichiro Nakamura.科学计算引论基于MATLAB的数值分析M.电子工业出版社,北京.2006.1,复习 分子物理学绘图,例4.6：利用气体分子运动的麦克斯韦速度分布率，求27C下氮分子运动的速度分布曲线，并求速度在300-500m/s范围内的分子所占的比例，讨论温度T及分子量对速度分布曲线的影响。 解： 建模 1.麦克斯韦速度分布率为： 2.考虑到

2、该公式较复杂，建立.m文件。,积分函数trapz(),程序： 1）子程序（mksw.m）： function f=mksw(T,mu,v) R=8.31; %气体常数 k=1.381*10(-23); %玻尔茨曼常数 NA=6.022*1023; %阿伏伽德罗数 m=mu/NA; %分子质量 f=4*pi*(m/(2*pi*k*T).(3/2) . .*exp(-m*v.2./(2*k*T).*V.2; %速度分布率,2）主程序：, T=input(绝对温度T=); mu=input(气体分子量mu=); vmin=input(速度下限vmin=); vmax=input(速度上限vmax=)

3、; v=0:1500; y=mksw(T,mu,v); plot(v,y); hold on; v1=vmin:vmax; %速度分布率 y1=mksw(T,mu,v1); fill(v1,500,300,y1,0,0,r); trapz(y1);,结果：,当输入: 绝对温度T=300 气体分子量mu=0.028 速度下限vmin=300 速度上限vmax=500 ans = 0.3763,复习我国古代数学家张丘在“算经”里提出一个世界数学史上有名的百鸡问题：鸡翁一，值钱五，鸡母一，值钱三，幼鸡三，值钱一，百钱买百鸡，问各几何？,解： 建模 (怎么建？),1）运用循环语句forend a.一重

4、； forifendend b.二重； for forifendend end c.三重。 for for forifendend end end 2）利用网格数组meshgrid( ),方法1-1：,for x=0:25 y=(100-7*x)/4; if mod(100-x-y,3)=0 end end,x+y+z=100 5x+3y+z/3=100 解方程得： 7x+4y=100 y=(100-7x)/4,百钱买百鸡,方法1-1：运用符号函数 syms+solve,syms x y z p=x;q=z; for y=0:33 f1=x+y+z; x=solve(f1-100); f2=5

5、*x+3*y+z/3; z=solve(f2-100); if mod(eval(z),3)=0 end,百钱买百鸡,方法1-2：,disp(鸡翁 鸡母 幼鸡); for x=0:20 for y=0:33 z=100-x-y; if 5*x+3*y+z/3=100 fprintf(%g %g %gn,x,y,z); end end end,鸡翁 鸡母 幼鸡 0 25 75 4 18 78 8 11 81 12 4 84,百钱买百鸡,方法1-3：,disp(鸡翁 鸡母 幼鸡); for x=0:20 for y=0:33 for z=3:3:99 if x+y+z=100 end end en

6、d end,百钱买百鸡,方法2-1：,x,y=meshgrid(0:20,0:33); t=(find(5*x+3*y+(100-x-y)/3=100); x(t) y(t) 100-x(t)-y(t),鸡翁： 0 4 8 12 鸡母： 25 18 11 4 幼鸡： 75 78 81 84,百钱买百鸡,一.单变量函数的计算和绘图 例5.1：已知 要求以0.01秒为间隔，绘出y及其导数的曲线. 分析：用diff(y,n)求Dy ， 每求导一次，y的维数减一。 Dy=diff(y)结果为Dy=y1-y2, 故 y=Dy/Dx= diff(y)/Dx,5.1 函数、极限和导数,b=0.1;t=0:b

7、:1.5;w= 4*sqrt(3); y=5*sqrt(3)*exp(-2*t).*sin(w*t+pi/3); plot(t,y); title(单变量绘图); xlabel(x);ylabel(y(t); grid on;hold on; Dy=diff(y)/b; plot(t(1:length(t)-1),Dy,*),程序：,plot(t(2:length(t),Dy,p) y1=-10*sqrt(3)*exp(-2*t) .*sin(w*t+pi/3)+60*exp(-2*t) .*cos(w*t+pi/3); hold on ; plot(t,y1,r) legend(y,Dy1,

8、 Dy2, y1),Dy和y1不重合呢？,例5.2.求点u=(1,2,3)到平面 3x-2y+z=4的距离,采用input()语句输入7个变量 烦 采用（1/2），为什么不用sqrt() 困 D输入4，为什么不是-4呢 晕,（提示：点(u,v,w)到面Ax+By+Cz+D=0公式为r=|Au+Bv+Cw+D|/(A2+B2+C2)1/2 ),u=input(点的坐标u=); v=input(直线方程系数v=); r=abs(u*v(1:3)+v(4)/sqrt(sum(v(1:3).2),点的坐标u=1,2,3 直线方程系数v=3,-2,1,-4 r = 0.5345,5.1 函数、极限和导数

9、,二.参变方程函数的计算和绘图 例5.3：已知炮弹初速为V0，发射角为alfan,画出其alfan为25，45，75度时的轨迹。 分析：,5.1 函数、极限和导数, V0=input(初始速度V0=); 初始速度V0=30 alfan=input(初始角alfan=); 初始角alfan=45 alfan1=alfan*pi/180; t=0:0.1:2*V0*sin(alfan1)/9.8; x=V0*cos(alfan1)*t; y=V0*sin(alfan1)*t-9.8*t.*t/2; plot(x,y); axis(equal),补充：三次抛物线方程y=ax3+cx，讨论参数a，c对

10、其图形的影响。,x=-3:.01:3; subplot(1,2,1); for c=-3:3 plot(x,x.3+c*x); hold on; end axis equal; axis(-3,3,-3,3),subplot(1,2,2); for a=-3:3 plot(x,a*x.3+x); hold on; end axis equal; axis(-3,3,-3,3),5.1 函数、极限和导数,三.极限 syms var1 var2 var3;多个符号变量函数 limit(f,x,a,s) 极限函数s=right,left或空 例5.4：求下列极限 1）, syms x; f=sin(

11、x)/x; limit(f,x,0),x=sym(x) 单个符号变量函数,5.1 函数、极限和导数, syms x; f=1/(x+2)-12/(x3+8) limit(f,x,-2) syms x a f=(a+x)3-a3)/x limit(f,x,0),2) 3),5.2 空间解析几何,一.求切点 例5.5:求曲线y=x3+3x-2上与直线y=100 x-1平行的切线的切点,并绘出曲线和切线。,建模：切点是其导数值为100的点 求导用函数diff() 求根用函数solve(),程序求切点 x=sym(x); y=x3+3*x-2; f=diff(y); x1=solve(f-100);

12、y1=x1.3+3*x1-2; c=y1-100*x1;hold on; plot(eval(x1),eval(y1),*), x=-100:0.1:100; y1=100*x+ eval(c(1); y2= 100*x+ eval(c(2); plot(x,y1,x,y2,r) y3=x.3+3*x-2; y4=100*x-1; plot(x,y3,x,y4,r) axis(-10,10,-500,500),程序绘图,5.2 空间解析几何,5.2 空间解析几何,二.求截面 例5.6:绘出用平行截面法后方程Z=x2-2y2构成的马鞍面形状。 建模：定义网格函数meshgrid() X,Y= m

13、eshgrid(x,y),5.2 空间解析几何,程序： clf x,y=meshgrid(-10:0.2:10); z1=(x.2-2*y.2); a=input(a=?); a=?20 z2=a*ones(size(x); subplot(1,2,1),mesh(x,y,z1); hold on;mesh(x,y,z2);,5.2 空间解析几何, v=-10,10,-10,10,-100,100; axis(v);grid; hold off; r0=abs(z1-z2) xx=r0.*x;yy=r0.*y;zz=r0.*z2; subplot(1,2,2), plot3(xx(r0=0),

14、yy(r0=0),zz(r0=0),*); axis(v);grid;,5.2 空间解析几何,5.2 空间解析几何,思考题： 绘出空间任意两曲面的交线。 （z1=x2-2y2 z2=2x-3y） 要求：方程为任意输入。,程序：,z1=input(输入方程z1=：,s); z2=input(输入方程z2=：,s); x,y=meshgrid(-10:0.2:10); z1=eval(z1); z2=eval(z2); mesh(x,y,z1); hold on; mesh(x,y,z2); r0=abs(z1-z2)=1; xx=r0.*x;yy=r0.*y;zz=r0.*z2; plot3(x

15、x(r0=0),yy(r0=0),zz(r0=0),*);,输入方程z1=：x.2-2*y.2 输入方程z2=：2*x-3*y,怎样改写源程序绘成下图？,5.3 数列和级数,例5.7:求数列 当x=1.2时的前6项。,程序： x=1.2; n=1:6; y=(-1).(n+1).*x.n./n 1.200 -0.720 0.576 -0.5184 0.498 -0.498,5.3 数列和级数,例5.8:求数列 当x=0.1,0.2,0.3,0.4时的前n项和与极限，并绘图说明。 程序： x=input(x=); x=0.1,0.2,0.3,0.4 n=input(n=); n=10, for

16、i=1:length(x) for j=1:n y(j,i)=(-1).(j+1).*x(i).j./j; if j1 y(j,i)=y(j-1,i)+y(j,i);end end end y,5.3 数列和级数,y =0.1000 0.2000 0.3000 0.4000 0.0950 0.1800 0.2550 0.3200 0.0953 0.1827 0.2640 0.3413 0.0953 0.1823 0.2620 0.3349 0.0953 0.1823 0.2625 0.3370 0.0953 0.1823 0.2623 0.3363 0.0953 0.1823 0.2624 0

17、.3365 0.0953 0.1823 0.2624 0.3365 0.0953 0.1823 0.2624 0.3365 0.0953 0.1823 0.2624 0.3365,5.3 数列和级数,syms n x=0.1,0.2,0.3,0.4; y=symsum(-1)(n+1)*x.n/n,1,inf) eval(y),symsum(a,n,n0,nn),y = log(11/10), log(6/5), log(13/10), log(7/5) ans = 0.0953 0.1823 0.2624 0.3365,5.3 数列和级数,x=0:0.01:5; y=symsum(-1)(n

18、+1)*x.n/n,1,inf); plot(x,eval(y),r)；plot(x,x,*),5.4 数值方法和符号积分,求解符号函数solve(eqn1,eqn2,.,var1,var2,.),一.求f(x)=0的解,5.4 数值方法和符号积分,例5.9:求曲线f(x)=x3-2x2-x+2的过零解，并绘图。 建模：求解符号函数是solve() 过零解:f(x)=0,5.4 数值方法和符号积分,程序： x=sym(x); y=x3-2*x2-x+2; z=solve(y); plot(eval(z),0,0,0,*r) hold on;x=-7:5; y=x.3-2*x.2-x+2; pl

19、ot(x,y),5.4 数值方法和符号积分,二.求定积分 int(f,x) int(f,x,a,b),5.4 数值方法和符号积分,例5.10:求椭球的体积。其方程为： 其中a,b,c为常数，用input（）输入。,5.4 数值方法和符号积分,程序： syms a b c z; f=pi*a*b*(c2-z2)/c2; V=int(f,z,-c,c) V =4/3*pi*a*b*c,建模： 用平面z=z0去截取椭球，其相交线是椭圆，其面积为,x,y,z=sphere(30); meshc(x*5,y*3,z*4),5.4 数值方法和符号积分,三.线性插值 interp1(x,y,xi,s),s：

20、linear/spline/cubic,程序：,x=0.0:0.25:1.0; y=0.9162,0.8109,0.6931,0.5596,0.4055; yi=0.9,0.7,0.6,0.5; xi=interp1(y,x,yi); yi,xi,5.4 数值方法和符号积分,四.曲线拟合 polyfit(x,y,n),x=0.1,0.4,0.5,0.7,0.7,0.9; y=0.61,0.92,0.99,1.52,1.47,2.03; c=polyfit(x,y,2); xx=x(1):0.1:x(length(x); yy=polyval(c,xx); plot(x,y,*,xx,yy),5

21、.5 线形代数,例5.11:用求逆矩阵的方法解线形方程组。 x+2y+3z=5 x+4y+9z=-2 x+8y+27z=6,建模： 1 2 3 5 x A= 1 4 9 b= -2 X= y 1 8 27 6 z,5.5 线形代数,即 AX=b 其解为X=A-1b 程序： A=1,2,3;1,4,9;1,8,27; b=5,-2,6; X=inv(a)*b %可换成：? 23.0000 -14.5000 3.6667,5.5 线形代数,例5.12:求矩阵的秩和迹，行列式，特征值和特征向量。 建模： 秩：行中最大线形无关组，rank(A) 迹：对角线元素之和，trace(A) 行列式：det(A

22、) 特征值和特征向量 :A=，, =eig(A),5.5 线形代数,程序： A=1,2,2;1,-1,1;4,-12,1; rank(A) ans = 3 trace(A) ans = 1 det(A) ans = 1,5.5 线形代数, v,d=eig(A) v = 0.9045 -0.7255 -0.7255 0.3015 -0.2176 - 0.0725i -0.2176 + 0.0725i -0.3015 0.5804 - 0.2902i 0.5804 + 0.2902i d = 1.0000 0 0 0 -0.0000 + 1.0000i 0 0 0 -0.0000 - 1.0000

23、i,1.求下列函数的极限：,2.求下列线形方程组的解，并求其系数矩阵的秩和迹，行列式，特征值和特征向量。 2x1+2x2- x3+ x4=4 4x1+3x2- x3+2x4=6 8x1+5x2-3x3+4x4=12 3x1+3x2-2x3+2x4=6 3.求x2+y2=z的体积(0=z=5)并绘立体图和z=3时的截面图。,练 习,郑州招聘会3万大学生拥挤图,谢 谢！,实验 部分答案,1.求任意两个正整数的最小公约数和最大公倍数。 解： 建模：,1）利用mod(x,y)循环求解。 2）利用mod(x,y)从大到小逐次求解。,程序1：,a=input(输入一个整数); b=input(输入另一个整

24、数); if ab m=a; a=b; b=m; end for i=b:-1:1 if rem(a,i)=0,程序2：,x=input(请输入任意两个正整数x=); y=max(x),min(x); while rem(y(1),y(2)=0 a=y(1); y(1)=y(2); y(2)=rem(a,y(2); end t=y(2),x(1)*x(2)/y(2); fprintf(任意两个正整数x=%f,%fn,x); fprintf(最大公约数和最小公倍数是%f %fn,t);,2. “同构数”是指这样的整数：它恰好出现在其平方数的右端。如：636，25625。试找出11000之间的全部

25、同构数。,解： 建模,1）运用循环语句forend 2）利用数组元素,方法1：,for n=1:1000 if rem(n2,10)=n|rem(n2,100)=n| rem(n2,1000)=n fprintf(%g ,n); end end 1 5 6 25 76 376 625,方法2：,for n=1:1000 i=1;m=n; while floor(m/10)=0 i=i+1; m=floor(m/10); end if rem(n2,10i)=n fprintf(%g ,n); end end 1 5 6 25 76 376 625,方法3：,x=1:999; y=rem(x(1:9).2,10),rem(x(10:99).2,100),rem(x(100:999).2,1000); t=(y=x); x(t) ans = 1 5 6 25 76 376 625,3. “完备数”是指一个数恰好等于它的因子之和。如6的因