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=); v=0
3、: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=500ans = 0.3763,复习我国古代数学家张丘在“算经”里提出一个世界数学史上有名的百鸡问题：鸡翁一，值钱五，鸡母一，值钱三，幼鸡三，值钱一，百钱买百鸡，问各几何？,解： 建模 (怎么建？),1）运用循环语句forend a.一重； fori
4、fendend b.二重； for forifendend end c.三重。 for for forifendend end end 2）利用网格数组meshgrid( ),方法1-1：,for x=0:25y=(100-7*x)/4;if mod(100-x-y,3)=0end 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:33f1=x+y+z;x=solve(f1-100);f2=5*x+3*y+z/3;z
5、=solve(f2-100);if mod(eval(z),3)=0 end,百钱买百鸡,方法1-2：,disp(鸡翁 鸡母 幼鸡); for x=0:20for y=0:33z=100-x-y;if 5*x+3*y+z/3=100fprintf(%g %g %gn,x,y,z);endend end,鸡翁 鸡母 幼鸡0 25 754 18 788 11 8112 4 84,百钱买百鸡,方法1-3：,disp(鸡翁 鸡母 幼鸡); for x=0:20for y=0:33for z=3:3:99if x+y+z=100endendend end,百钱买百鸡,方法2-1：,x,y=meshgri
6、d(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:1.5;w= 4*sqrt(3); y=5*sqrt(3)*exp
7、(-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, Dy2, y1),Dy和y1不重合呢？,例5.2.求点u=(1,2,3
8、)到平面 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 函数、极限和导数,二.参变方程函数的计算和绘图例5.3：已知炮弹初速为V0，发射角为alfan
9、,画出其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对其图形的影响。,x=-3:.01:3; subplot(1,2,1); for c
10、=-3:3plot(x,x.3+c*x);hold on; end axis equal; axis(-3,3,-3,3),subplot(1,2,2); for a=-3:3plot(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(x)/x; limit(f,x,0),x=sym(x) 单个符号变量函数,5.1 函数、极限和
11、导数, 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=100x-1平行的切线的切点,并绘出曲线和切线。,建模：切点是其导数值为100的点求导用函数diff()求根用函数solve(),程序求切点 x=sym(x); y=x3+3*x-2; f=diff(y); x1=solve(f-100); y1=x1.3+3*x1-2; c=y1-100*x1;hold on; plot(eval(x1),
12、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= meshgrid(x,y),5.2 空间解析几何,程序： clf x,y=meshgrid(-10:0.2:
13、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),yy(r0=0),zz(r0=0),*); axis(v);grid;,5.2 空间解析几何,5.2 空间
14、解析几何,思考题： 绘出空间任意两曲面的交线。 （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(xx(r0=0),yy(r0=0),zz(r0=0),*);,怎样改写源程序绘成下图？,5.3 数列和级数,例
15、5.7:求数列 当x=1.2时的前6项。,程序： x=1.2; n=1:6; y=(-1).(n+1).*x.n./n1.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 i=1:length(x)for j=1:ny(j,i)=(-1).(j+1).*x(i).j./j;if j1 y(j,i)=y(j-1,i)+y(j,i);end endend y
16、,5.3 数列和级数,y =0.1000 0.2000 0.3000 0.40000.0950 0.1800 0.2550 0.32000.0953 0.1827 0.2640 0.34130.0953 0.1823 0.2620 0.33490.0953 0.1823 0.2625 0.33700.0953 0.1823 0.2623 0.33630.0953 0.1823 0.2624 0.33650.0953 0.1823 0.2624 0.33650.0953 0.1823 0.2624 0.33650.0953 0.1823 0.2624 0.3365,5.3 数列和级数,syms
17、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+1)*x.n/n,1,inf); plot(x,eval(y),r)；plot(x,x,*),5.4 数值方法和符号积分,求解符号函数solve(eqn1,eqn2,.,var1,var2,.),一.求f(x
18、)=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; plot(x,y),5.4 数值方法和符号积分,二.求定积分int(f,x) int(f,x,a,b),5.4 数值方法和符号积分,例5.10:求椭球的体积。其方程为：其中a,b,c为常数，用input（）输入。,5.
19、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：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
20、,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.5 线形代数,例5.11:用求逆矩阵的方法解线形方程组。x+2y+3z=5x+4y+9z=-2x+8y+27z=6,建模：1 2 3 5 x A= 1 4 9 b= -2 X= y1 8 27 6 z,5.5 线形代数
21、,即 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) 特征值和特征向量 :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 线形代数
22、, 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.0000i,1.求下列函数的极限：,2.求下列线形方程组的解，并求其系数矩阵的秩和迹，行列式，特征值和特征向量。2x1+2x2- x3+ x4=44x1+3x2- x3+2x4=68x1+5x2-3x3+4x4=123x1+3x2-2x3+2x4=6
23、 3.求x2+y2=z的体积(0=z=5)并绘立体图和z=3时的截面图。,练 习,郑州招聘会3万大学生拥挤图,谢 谢！,实验 部分答案,1.求任意两个正整数的最小公约数和最大公倍数。解：建模：,1）利用mod(x,y)循环求解。 2）利用mod(x,y)从大到小逐次求解。,程序1：,a=input(输入一个整数); b=input(输入另一个整数); if abm=a; a=b; b=m; end for i=b:-1:1if rem(a,i)=0,程序2：,x=input(请输入任意两个正整数x=); y=max(x),min(x); while rem(y(1),y(2)=0a=y(1);
24、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之间的全部同构数。,解： 建模,1）运用循环语句forend2）利用数组元素,方法1：,for n=1:1000if rem(n2,10)=n|rem(n2,100)=n|rem(n2,1000)=nfprintf(%g ,n);end end1 5 6 25 76 376
25、 625,方法2：,for n=1:1000i=1;m=n;while floor(m/10)=0i=i+1;m=floor(m/10);endif rem(n2,10i)=nfprintf(%g ,n);end end1 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的因子为1，2，3，而6123。编制程序找出1500之间
26、的全部完备数。,解： 建模1）运用累加计数:y=y+2）利用数组元素添加：y=y,方法1：,for i=1:500y=0;for j=1:i/2if mod(i,j)=0y=y+j;endendif i=yfprintf(%d ,i);end end6 28 496,方法2：,for i=1:500y=;for j=1:i/2if mod(i,j)=0y=y,j;endendif i=sum(y)fprintf(%d ,i);end end6 28 496,4. “自然数对”是指两个自然数的和与差都是平方数。如8与17的和81725与其差1789都是平方数，则8和17称自然数对。求出1100间
27、的自然数对。,解： 建模关键是如何判断平方数： m=fix(m),程序：,for i=1:100for j=i+1:100m=sqrt(j+i);n=sqrt(j-i);if m=fix(m)endend end,5.求x2+y2=z的体积(0=z=5)并绘立体图和z=3时的截面图。,syms z; f=pi*z; V=int(f,z,0,5)V =25/2*pi,用平面z=z0去截取该立体，其相交线是圆，其面积为：S= (x2+y2)1/2)2= z,?该立体不是圆锥体,故不能用公式V= r2h/3,求体积：,画立体图和横截面,clear all x,y=meshgrid(-3:0.2:3); z1=(x.2+y.2); mesh(x,y,z1) z2=3*ones(size(x); hold on; mesh(x,y,z2),（-3:0.2:3 ，-3:0.2:3）length(x),区别？,画截线： r0=abs(z1-z2)=1; xx=r0.*x;yy=r0.*y;zz=r0.*z2; plot3(xx(r0=0),yy(r0=0),zz(r0=0),*);,
