MATLAB程序设计与应用-5

第五讲MATLAB数值计算符号计算授课教师：田鹏数值计算+符号计算5.1函数、极限和导数5.2空间解析几何5.3数列和级数5.4数值方法和符号积分5.5线形代数参考 书 关于书的成语关于读书的排比句社区图书漂流公约怎么写关于读书的小报汉书pdf 目：ShoichiroNakamura.科学计算引论——基于MATLAB的数值 分析 定性数据统计分析pdf销售业绩分析模板建筑结构震害分析销售进度分析表京东商城竞争战略分析 [M].电子工业出版社,北京.2006.1复习分子物理学绘图例4.6：利用气体分子运动的麦克斯韦速度分布率，求27°C下氮分子运动的速度分布曲线，并求速度在300-500m/s范围内的分子所占的比例，讨论温度T及分子量对速度分布曲线的影响。解：建模1.麦克斯韦速度分布率为：2.考虑到该公式较复杂，建立.m文件。积分函数trapz()程序：1）子程序（mksw.m）：functionf=mksw(T,mu,v)R=8.31;%气体常数k=1.381*10^(-23);%玻尔茨曼常数NA=6.022*10^23;%阿伏伽德罗数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:1500;>>y=mksw(T,mu,v);>>plot(v,y);holdon;>>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）运用循环语句for…enda.一重；for…if…end…endb.二重；for…for…if…end…end…endc.三重。for…for…for…if…end…end…end…end2）利用网格数组meshgrid()方法1-1：forx=0:25y=(100-7*x)/4;ifmod(100-x-y,3)==0&y>=0fprintf('%g%g%g\n',x,y,100-x-y);endendx+y+z=1005x+3y+z/3=100解方程得：7x+4y=100y=(100-7x)/4百钱买百鸡方法1-1：运用符号函数syms+solvesymsxyzp=x;q=z;fory=0:33f1=x+y+z;x=solve(f1-100);f2=5*x+3*y+z/3;z=solve(f2-100);ifmod(eval(z),3)==0&y+eval(z)<=100fprintf('%g%g%g\n',100-y-eval(z),y,eval(z));endx=p;z=q;end为什么不用z=eval(solve(f2-100));百钱买百鸡方法1-2：disp('鸡翁鸡母幼鸡');forx=0:20fory=0:33z=100-x-y;if5*x+3*y+z/3==100fprintf('%g%g%g\n',x,y,z);endendend鸡翁鸡母幼鸡02575418788118112484百钱买百鸡方法1-3：disp('鸡翁鸡母幼鸡');forx=0:20fory=0:33forz=3:3:99ifx+y+z==100&5*x+3*y+z/3==100fprintf('%g%g%g\n',x,y,z);endendendend百钱买百鸡方法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)'鸡翁：04812鸡母：2518114幼鸡：75788184百钱买百鸡一.单变量函数的计算和绘图例5.1：已知要求以0.01秒为间隔，绘出y及其导数的曲线.分析：用diff(y,n)求Dy，每求导一次，y的维数减一。Dy=diff(y)结果为Dy=y1-y2,故y’=Dy/Dx=diff(y)/Dx5.1函数、极限和导数b=0.1;t=0:b: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)');gridon;holdon;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);holdon;plot(t,y1,‘r')legend('y','Dy1','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.53455.1函数、极限和导数二.参变方程函数的计算和绘图例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对其图形的影响。x=-3:.01:3;subplot(1,2,1);forc=-3:3plot(x,x.^3+c*x);holdon;endaxisequal;axis([-3,3,-3,3])subplot(1,2,2);fora=-3:3plot(x,a*x.^3+x);holdon;endaxisequal;axis([-3,3,-3,3])5.1函数、极限和导数三.极限symsvar1var2var3……;多个符号变量函数limit(f,x,a,’s’)极限函数s=right,left或空例5.4：求下列极限1）>>symsx;>>f=sin(x)/x;>>limit(f,x,0)x=sym(‘x’)单个符号变量函数5.1函数、极限和导数>>symsx;>>f=1/(x+2)-12/(x^3+8)>>limit(f,x,-2)>>symsxa>>f=((a+x)^3-a^3)/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=x^3+3*x-2;>>f=diff(y);>>x1=solve(f-100);>>y1=x1.^3+3*x1-2;>>c=y1-100*x1;holdon;>>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])程序——绘图y1=100*x+eval(c);plot(x,y1,'r')5.2空间解析几何5.2空间解析几何二.求截面例5.6:绘出用平行截面法后方程Z=x2-2y2构成的马鞍面形状。建模：定义网格函数meshgrid()[X,Y]=meshgrid(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);>>holdon;mesh(x,y,z2);5.2空间解析几何>>v=[-10,10,-10,10,-100,100];>>axis(v);grid;holdoff;>>r0=abs(z1-z2)<=1;>>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空间解析几何思考题：绘出空间任意两曲面的交线。（z1=x2-2y2z2=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);holdon;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),'*');输入方程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./n1.200-0.7200.576-0.51840.498-0.4985.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>>fori=1:length(x)forj=1:ny(j,i)=(-1).^(j+1).*x(i).^j./j;ifj>1y(j,i)=y(j-1,i)+y(j,i);endendend>>y5.3数列和级数y=0.10000.20000.30000.40000.09500.18000.25500.32000.09530.18270.26400.34130.09530.18230.26200.33490.09530.18230.26250.33700.09530.18230.26230.33630.09530.18230.26240.33650.09530.18230.26240.33650.09530.18230.26240.33650.09530.18230.26240.33655.3数列和级数symsnx=[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.09530.18230.26240.33655.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)=0的解例X2+xy+y=3X2-4y+3=0[x,y]=solve('x^2+x*y+y=3',…'x^2-4*y+3=0')5.4数值方法和符号积分例5.9:求曲线f(x)=x3-2x2-x+2的过零解，并绘图。建模：求解符号函数是solve()过零解:f(x)=05.4数值方法和符号积分程序：>>x=sym('x');>>y=x^3-2*x^2-x+2;>>z=solve(y);>>plot(eval(z),[0,0,0],'*r')>>holdon;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.4数值方法和符号积分程序：>>symsabcz;>>f=pi*a*b*(c^2-z^2)/c^2;>>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’)0.40551.00.50.55960.750.60.69310.50.70.81090.250.90.91620.0yiyxs：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)2.030.91.470.71.520.70.990.50.920.40.610.1yxx=[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建模：1235xA=149b=-2X=y18276z5.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.50003.66675.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=15.5线形代数>>[v,d]=eig(A)v=0.9045-0.7255-0.72550.3015-0.2176-0.0725i-0.2176+0.0725i-0.30150.5804-0.2902i0.5804+0.2902id=1.0000000-0.0000+1.0000i000-0.0000-1.0000i1.求下列函数的极限：2.求下列线形方程组的解，并求其系数矩阵的秩和迹，行列式，特征值和特征向量。2x1+2x2-x3+x4=44x1+3x2-x3+2x4=68x1+5x2-3x3+4x4=123x1+3x2-2x3+2x4=63.求x2+y2=z的体积(0<=z<=5)并绘立体图和z=3时的截面图。练习郑州招聘会3万大学生拥挤图谢谢！实验部分 答案 八年级地理上册填图题岩土工程勘察试题省略号的作用及举例应急救援安全知识车间5s试题及答案 1.求任意两个正整数的最小公约数和最大公倍数。解：建模：1）利用mod(x,y)循环求解。2）利用mod(x,y)从大到小逐次求解。程序1：a=input('输入一个整数');b=input('输入另一个整数');ifa