MATLAB数学实验课后答案

Fourshortwordssumupwhathasliftedmostsuccessfulindividualsabovethecrowd:alittlebitmore.------------------------------------------author------------------------------------------dateMATLAB数学实验课后答案MATLAB数学实验课后答案MATLAB数学实验课后答案----------------------------------------------------------------------------------------------------MATLAB数学实验课后答案--------------------------------------------------数学实验MATLAB参考答案（重要部分）P20,ex1(5)等于[exp(1),exp(2);exp(3),exp(4)](7)3=1*3,8=2*4(8)a为各列最小值，b为最小值所在的行号(10)1>=4,false,2>=3,false,3>=2,ture,4>=1,ture(11)答案表明：编址第2元素满足不等式(30>=20)和编址第4元素满足不等式(40>=10)(12)答案表明：编址第2行第1列元素满足不等式(30>=20)和编址第2行第2列元素满足不等式(40>=10)P20,ex2(1)a,b,c的值尽管都是1，但数据类型分别为数值，字符，逻辑，注意a与c相等，但他们不等于b(2)double(fun)输出的分别是字符a,b,s,(,x,)的ASCII码P20,ex3>>r=2;p=0.5;n=12;>>T=log(r)/n/log(1+0.01*p)T=11.5813P20,ex4>>x=-2:0.05:2;f=x.^4-2.^x;>>[fmin,min_index]=min(f)fmin=-1.3907%最小值min_index=54%最小值点编址>>x(min_index)ans=0.6500%最小值点>>[f1,x1_index]=min(abs(f))%求近似根--绝对值最小的点f1=0.0328x1_index=24>>x(x1_index)ans=-0.8500>>x(x1_index)=[];f=x.^4-2.^x;%删去绝对值最小的点以求函数绝对值次小的点>>[f2,x2_index]=min(abs(f))%求另一近似根--函数绝对值次小的点f2=0.0630x2_index=65>>x(x2_index)ans=1.2500P20,ex5>>z=magic(10)z=929918156774515840988071416735557644148188202254566370478587192136062697128869325296168755234172476839042492633652358289914830323966796139597293138457210129496783537444653111810077843643502759>>sum(z)ans=505505505505505505505505505505>>sum(diag(z))ans=505>>z(:,2)/sqrt(3)ans=57.157746.188046.765450.229553.693613.85642.88683.46416.928210.3923>>z(8,:)=z(8,:)+z(3,:)z=9299181567745158409880714167355576441481882022545663704785871921360626971288693252961687552341724768390424926336523582899148303239668387101115119838710111511910129496783537444653111810077843643502759P40ex1先在编辑器窗口写下列M函数，保存为eg2_1.mfunction[xbar,s]=ex2_1(x)n=length(x);xbar=sum(x)/n;s=sqrt((sum(x.^2)-n*xbar^2)/(n-1));例如>>x=[81706551766690876177];>>[xbar,s]=ex2_1(x)xbar=72.4000s=12.1124P40ex2s=log(1);n=0;whiles<=100n=n+1;s=s+log(1+n);endm=n计算结果m=37P40ex3clear;F(1)=1;F(2)=1;k=2;x=0;e=1e-8;a=(1+sqrt(5))/2;whileabs(x-a)>ek=k+1;F(k)=F(k-1)+F(k-2);x=F(k)/F(k-1);enda,x,k计算至k=21可满足精度P40ex4clear;tic;s=0;fori=1:1000000s=s+sqrt(3)/2^i;ends,toctic;s=0;i=1;whilei<=1000000s=s+sqrt(3)/2^i;i=i+1;ends,toctic;s=0;i=1:1000000;s=sqrt(3)*sum(1./2.^i);s,tocP40ex5t=0:24;c=[15141414141516182022232528...313231292725242220181716];plot(t,c)P40ex6(1)clear;fplot('x^2*sin(x^2-x-2)',[-2,2])x=-2:0.1:2;y=x.^2.*sin(x.^2-x-2);plot(x,y)y=inline('x^2*sin(x^2-x-2)');fplot(y,[-22])(2)参数方法t=linspace(0,2*pi,100);x=2*cos(t);y=3*sin(t);plot(x,y)(3)x=-3:0.1:3;y=x;[x,y]=meshgrid(x,y);z=x.^2+y.^2;surf(x,y,z)(4)x=-3:0.1:3;y=-3:0.1:13;[x,y]=meshgrid(x,y);z=x.^4+3*x.^2+y.^2-2*x-2*y-2*x.^2.*y+6;surf(x,y,z)(5)t=0:0.01:2*pi;x=sin(t);y=cos(t);z=cos(2*t);plot3(x,y,z)(6)theta=linspace(0,2*pi,50);fai=linspace(0,pi/2,20);[theta,fai]=meshgrid(theta,fai);x=2*sin(fai).*cos(theta);y=2*sin(fai).*sin(theta);z=2*cos(fai);surf(x,y,z)(7)x=linspace(0,pi,100);y1=sin(x);y2=sin(x).*sin(10*x);y3=-sin(x);plot(x,y1,x,y2,x,y3)page41,ex7x=-1.5:0.05:1.5;y=1.1*(x>1.1)+x.*(x<=1.1).*(x>=-1.1)-1.1*(x<-1.1);plot(x,y)page41,ex8分别使用whichtrapz,typetrapz,dirC:\MATLAB7\toolbox\matlab\datafun\page41,ex9clear;close;x=-2:0.1:2;y=x;[x,y]=meshgrid(x,y);a=0.5457;b=0.7575;p=a*exp(-0.75*y.^2-3.75*x.^2-1.5*x).*(x+y>1);p=p+b*exp(-y.^2-6*x.^2).*(x+y>-1).*(x+y<=1);p=p+a*exp(-0.75*y.^2-3.75*x.^2+1.5*x).*(x+y<=-1);mesh(x,y,p)page41,ex10lookforlyapunovhelplyap>>A=[123;456;780];C=[2-5-22;-5-24-56;-22-56-16];>>X=lyap(A,C)X=1.0000-1.0000-0.0000-1.00002.00001.0000-0.00001.00007.0000Chapter3%Exercise1>>a=[1,2,3];b=[2,4,3];a./b,a.\b,a/b,a\bans=0.50000.50001.0000ans=221ans=0.6552%一元方程组x[2,4,3]=[1,2,3]的近似解ans=0000000.66671.33331.0000%矩阵方程[1,2,3][x11,x12,x13;x21,x22,x23;x31,x32,x33]=[2,4,3]的特解Exercise2(1)>>A=[41-1;32-6;1-53];b=[9;-2;1];>>rank(A),rank([A,b])%[A,b]为增广矩阵ans=3ans=3%可见方程组唯一解>>x=A\bx=2.38301.48942.0213Exercise2(2)>>A=[4-33;32-6;1-53];b=[-1;-2;1];>>rank(A),rank([A,b])ans=3ans=3%可见方程组唯一解>>x=A\bx=-0.4706-0.29410Exercise2(3)>>A=[41;32;1-5];b=[1;1;1];>>rank(A),rank([A,b])ans=2ans=3%可见方程组无解>>x=A\bx=0.3311-0.1219%最小二乘近似解Exercise2(4)>>a=[2,1,-1,1;1,2,1,-1;1,1,2,1];b=[123]';%注意b的写法>>rank(a),rank([a,b])ans=3ans=3%rank(a)==rank([a,b])<4说明有无穷多解>>a\bans=1010%一个特解Exercise3>>a=[2,1,-1,1;1,2,1,-1;1,1,2,1];b=[1,2,3]';>>x=null(a),x0=a\bx=-0.62550.6255-0.20850.4170x0=1010%通解kx+x0Exercise4>>x0=[0.20.8]';a=[0.990.05;0.010.95];>>x1=a*x,x2=a^2*x,x10=a^10*x>>x=x0;fori=1:1000,x=a*x;end,xx=0.83330.1667>>x0=[0.80.2]';>>x=x0;fori=1:1000,x=a*x;end,xx=0.83330.1667>>[v,e]=eig(a)v=0.9806-0.70710.19610.7071e=1.0000000.9400>>v(:,1)./xans=1.17671.1767%成比例，说明x是最大特征值对应的特征向量Exercise5%用到公式(3.11)(3.12)>>B=[6,2,1;2.25,1,0.2;3,0.2,1.8];x=[25520]';>>C=B/diag(x)C=0.24000.40000.05000.09000.20000.01000.12000.04000.0900>>A=eye(3,3)-CA=0.7600-0.4000-0.0500-0.09000.8000-0.0100-0.1200-0.04000.9100>>D=[171717]';x=A\Dx=37.569625.786224.7690%Exercise6(1)>>a=[41-1;32-6;1-53];det(a),inv(a),[v,d]=eig(a)ans=-94ans=0.2553-0.02130.04260.1596-0.1383-0.22340.1809-0.2234-0.0532v=0.0185-0.9009-0.3066-0.7693-0.1240-0.7248-0.6386-0.41580.6170d=-3.05270003.67600008.3766%Exercise6(2)>>a=[11-1;02-1;-120];det(a),inv(a),[v,d]=eig(a)ans=1ans=2.0000-2.00001.00001.0000-1.00001.00002.0000-3.00002.0000v=-0.57730.5774+0.0000i0.5774-0.0000i-0.57730.57740.5774-0.57740.5773-0.0000i0.5773+0.0000id=1.00000001.0000+0.0000i0001.0000-0.0000i%Exercise6(3)>>A=[5765;71087;68109;57910]A=5765710876810957910>>det(A),inv(A),[v,d]=eig(A)ans=1ans=68.0000-41.0000-17.000010.0000-41.000025.000010.0000-6.0000-17.000010.00005.0000-3.000010.0000-6.0000-3.00002.0000v=0.83040.09330.39630.3803-0.5016-0.30170.61490.5286-0.20860.7603-0.27160.55200.1237-0.5676-0.62540.5209d=0.010200000.843100003.8581000030.2887%Exercise6(4)、(以n=5为例)%关键是矩阵的定义%方法一（三个for）n=5;fori=1:n,a(i,i)=5;endfori=1:(n-1),a(i,i+1)=6;endfori=1:(n-1),a(i+1,i)=1;enda%方法二（一个for）n=5;a=zeros(n,n);a(1,1:2)=[56];fori=2:(n-1),a(i,[i-1,i,i+1])=[156];enda(n,[n-1n])=[15];a%方法三（不用for）n=5;a=diag(5*ones(n,1));b=diag(6*ones(n-1,1));c=diag(ones(n-1,1));a=a+[zeros(n-1,1),b;zeros(1,n)]+[zeros(1,n);c,zeros(n-1,1)]%下列计算>>det(a)ans=665>>inv(a)ans=0.3173-0.58651.0286-1.62411.9489-0.09770.4887-0.85711.3534-1.62410.0286-0.14290.5429-0.85711.0286-0.00750.0376-0.14290.4887-0.58650.0015-0.00750.0286-0.09770.3173>>[v,d]=eig(a)v=-0.7843-0.7843-0.92370.9860-0.92370.5546-0.5546-0.3771-0.00000.3771-0.2614-0.26140.0000-0.16430.00000.0924-0.09240.0628-0.0000-0.0628-0.0218-0.02180.02570.02740.0257d=0.7574000009.2426000007.4495000005.0000000002.5505%Exercise7(1)>>a=[41-1;32-6;1-53];[v,d]=eig(a)v=0.0185-0.9009-0.3066-0.7693-0.1240-0.7248-0.6386-0.41580.6170d=-3.05270003.67600008.3766>>det(v)ans=-0.9255%v行列式正常,特征向量线性相关，可对角化>>inv(v)*a*v%验算ans=-3.05270.0000-0.00000.00003.6760-0.0000-0.0000-0.00008.3766>>[v2,d2]=jordan(a)%也可用jordanv2=0.07980.00760.91270.1886-0.31410.1256-0.1605-0.26070.4213%特征向量不同d2=8.3766000-3.0527-0.0000i0003.6760+0.0000i>>v2\a*v2ans=8.376600.00000.0000-3.05270.00000.00000.00003.6760>>v(:,1)./v2(:,2)%对应相同特征值的特征向量成比例ans=2.44912.44912.4491%Exercise7(2)>>a=[11-1;02-1;-120];[v,d]=eig(a)v=-0.57730.5774+0.0000i0.5774-0.0000i-0.57730.57740.5774-0.57740.5773-0.0000i0.5773+0.0000id=1.00000001.0000+0.0000i0001.0000-0.0000i>>det(v)ans=-5.0566e-028-5.1918e-017i%v的行列式接近0,特征向量线性相关，不可对角化>>[v,d]=jordan(a)v=1011001-10d=110011001%jordan标准形不是对角的，所以不可对角化%Exercise7(3)>>A=[5765;71087;68109;57910]A=5765710876810957910>>[v,d]=eig(A)v=0.83040.09330.39630.3803-0.5016-0.30170.61490.5286-0.20860.7603-0.27160.55200.1237-0.5676-0.62540.5209d=0.010200000.843100003.8581000030.2887>>inv(v)*A*vans=0.01020.0000-0.00000.00000.00000.8431-0.0000-0.0000-0.00000.00003.8581-0.0000-0.0000-0.0000030.2887%本题用jordan不行,原因未知%Exercise7(4)参考6(4)和7(1),略%Exercise8只有(3)对称,且特征值全部大于零,所以是正定矩阵.%Exercise9(1)>>a=[4-313;2-135;1-1-1-1;3-234;7-6-70]>>rank(a)ans=3>>rank(a(1:3,:))ans=2>>rank(a([124],:))%1,2,4行为最大无关组ans=3>>b=a([124],:)';c=a([35],:)';>>b\c%线性表示的系数ans=0.50005.0000-0.50001.00000-5.0000%Exercise10>>a=[1-22;-2-24;24-2]>>[v,d]=eig(a)v=0.33330.9339-0.12930.6667-0.3304-0.6681-0.66670.1365-0.7327d=-7.00000002.00000002.0000>>v'*vans=1.00000.00000.00000.00001.000000.000001.0000%v确实是正交矩阵%Exercise11%设经过6个电阻的电流分别为i1,...,i6.列方程组如下%20-2i1=a;5-3i2=c;a-3i3=c;a-4i4=b;c-5i5=b;b-3i6=0;%i1=i3+i4;i5=i2+i3;i6=i4+i5;%计算如下>>A=[100200000;001030000;10-100-3000;1-10000-400;0-110000-50;01000000-3;00010-1-100;0000-1-1010;000000-1-11];>>b=[2050000000]';A\bans=13.34536.44018.54203.3274-1.18071.60111.72630.42042.1467%Exercise12>>A=[123;456;780];>>left=sum(eig(A)),right=sum(trace(A))left=6.0000right=6>>left=prod(eig(A)),right=det(A)%原题有错,(-1)^n应删去left=27.0000right=27>>fA=(A-p(1)*eye(3,3))*(A-p(2)*eye(3,3))*(A-p(3)*eye(3,3))fA=1.0e-012*0.08530.14210.02840.14210.14210-0.0568-0.11370.1705>>norm(fA)%f(A)范数接近0ans=2.9536e-013%Exercise1(1)roots([111])%Exercise1(2)roots([30-402-1])%Exercise1(3)p=zeros(1,24);p([1171822])=[5-68-5];roots(p)%Exercise1(4)p1=[23];p2=conv(p1,p1);p3=conv(p1,p2);p3(end)=p3(end)-4;%原p3最后一个分量-4roots(p3)%Exercise2fun=inline('x*log(sqrt(x^2-1)+x)-sqrt(x^2-1)-0.5*x');fzero(fun,2)】%Exercise3fun=inline('x^4-2^x');fplot(fun,[-22]);gridon;fzero(fun,-1),fzero(fun,1),fminbnd(fun,0.5,1.5)%Exercise4fun=inline('x*sin(1/x)','x');fplot(fun,[-0.10.1]);x=zeros(1,10);fori=1:10,x(i)=fzero(fun,(i-0.5)*0.01);end;x=[x,-x]%Exercise5fun=inline('[9*x(1)^2+36*x(2)^2+4*x(3)^2-36;x(1)^2-2*x(2)^2-20*x(3);16*x(1)-x(1)^3-2*x(2)^2-16*x(3)^2]','x');[a,b,c]=fsolve(fun,[000])%Exercise6fun=@(x)[x(1)-0.7*sin(x(1))-0.2*cos(x(2)),x(2)-0.7*cos(x(1))+0.2*sin(x(2))];[a,b,c]=fsolve(fun,[0.50.5])%Exercise7clear;close;t=0:pi/100:2*pi;x1=2+sqrt(5)*cos(t);y1=3-2*x1+sqrt(5)*sin(t);x2=3+sqrt(2)*cos(t);y2=6*sin(t);plot(x1,y1,x2,y2);gridon;%作图发现4个解的大致位置，然后分别求解y1=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[1.5,2])y2=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[1.8,-2])y3=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[3.5,-5])y4=fsolve('[(x(1)-2)^2+(x(2)-3+2*x(1))^2-5,2*(x(1)-3)^2+(x(2)/3)^2-4]',[4,-4])%Exercise8(1)clear;fun=inline('x.^2.*sin(x.^2-x-2)');fplot(fun,[-22]);gridon;%作图观察x(1)=-2;x(3)=fminbnd(fun,-1,-0.5);x(5)=fminbnd(fun,1,2);fun2=inline('-x.^2.*sin(x.^2-x-2)');x(2)=fminbnd(fun2,-2,-1);x(4)=fminbnd(fun2,-0.5,0.5);x(6)=2feval(fun,x)%答案:以上x(1)(3)(5)是局部极小，x(2)(4)(6)是局部极大，从最后一句知道x(1)全局最小，x(2)最大。%Exercise8(2)clear;fun=inline('3*x.^5-20*x.^3+10');fplot(fun,[-33]);gridon;%作图观察x(1)=-3;x(3)=fminsearch(fun,2.5);fun2=inline('-(3*x.^5-20*x.^3+10)');x(2)=fminsearch(fun2,-2.5);x(4)=3;feval(fun,x)%Exercise8(3)fun=inline('abs(x^3-x^2-x-2)');fplot(fun,[03]);gridon;%作图观察fminbnd(fun,1.5,2.5)fun2=inline('-abs(x^3-x^2-x-2)');fminbnd(fun2,0.5,1.5)%Exercise9close;x=-2:0.1:1;y=-7:0.1:1;[x,y]=meshgrid(x,y);z=y.^3/9+3*x.^2.*y+9*x.^2+y.^2+x.*y+9;mesh(x,y,z);gridon;%作图观察fun=inline('x(2)^3/9+3*x(1)^2*x(2)+9*x(1)^2+x(2)^2+x(1)*x(2)+9');x=fminsearch(fun,[00])%求极小值fun2=inline('-(x(2)^3/9+3*x(1)^2*x(2)+9*x(1)^2+x(2)^2+x(1)*x(2)+9)');x=fminsearch(fun2,[0-5])%求极大值%Exercise10clear;t=0:24;c=[15141414141516182022232528...313231292725242220181716];p2=polyfit(t,c,2)p3=polyfit(t,c,3)fun=inline('a(1)*exp(a(2)*(t-14).^2)','a','t');a=lsqcurvefit(fun,[00],t,c)%初值可以试探f=feval(fun,a,t)norm(f-c)%拟合效果plot(t,c,t,f)%作图检验fun2=inline('b(1)*sin(pi/12*t+b(2))+20','b','t');%原题修改f(x)+20b=lsqcurvefit(fun2,[00],t,c)figuref2=feval(fun2,b,t)norm(f2-c)%拟合效果plot(t,c,t,f2)%作图检验%Exercise11fun=inline('(1-x)*sqrt(10.52+x)-3.06*x*sqrt(1+x)*sqrt(5)');x=fzero(fun,0,1)%Exercise12r=5.04/12/100;N=20*12;x=7500*180%房屋总价格y=x*0.3%首付款额x0=x-y%贷款总额a=(1+r)^N*r*x0/((1+r)^N-1)%月付还款额r1=4.05/12/100;x1=10*10000;%公积金贷款a1=(1+r1)^N*r1*x1/((1+r1)^N-1)x2=x0-x1%商业贷款a2=(1+r)^N*r*x2/((1+r)^N-1)a=a1+a2%Exercise13%列方程th*R^2+(pi-2*th)*r^2-R*r*sin(th)=pi*r^2/2%化简得sin(2*th)-2*th*cos(2*th)=pi/2%以下Matlab计算clear;fun=inline('sin(2*th)-2*th*cos(2*th)-pi/2','th')th=fsolve(fun,pi/4)R=20*cos(th)%Exercise14%先在Editor窗口写M函数保存functionx=secant(fname,x0,x1,e)whileabs(x0-x1)>e,x=x1-(x1-x0)*feval(fname,x1)/(feval(fname,x1)-feval(fname,x0));x0=x1;x1=x;end%再在指令窗口fun=inline('x*log(sqrt(x^2-1)+x)-sqrt(x^2-1)-0.5*x');secant(fun,1,2,1e-8)%Exercise15%作系数为a,初值为xo,从第m步到第n步迭代过程的M函数：functionf=ex4_15fun(a,x0,m,n)x(1)=x0;y(1)=a*x(1)+1;x(2)=y(1);ifm<2,plot([x(1),x(1),x(2)],[0,y(1),y(1)]);holdon;endfori=2:ny(i)=a*x(i)+1;x(i+1)=y(i);ifi>m,plot([x(i),x(i),x(i+1)],[y(i-1),y(i),y(i)]);endendholdoff;%M脚本文件subplot(2,2,1);ex4_15fun(0.9,1,1,20);subplot(2,2,2);ex4_15fun(-0.9,1,1,20);subplot(2,2,3);ex4_15fun(1.1,1,1,20);subplot(2,2,4);ex4_15fun(-1.1,1,1,20);%Exercise16%设夹角t,问题转化为minf=5/sin(t)+10/cos(t)%取初始值pi/4,计算如下fun=@(t)5/sin(t)+10/cos(t);[t,f]=fminsearch(fun,pi/4)t=0.6709f=20.8097%Exercise17%提示：x(k+2)=f(x(k))=a^2*x(k)*(1-x(k))*(1-a*x(k)*(1-x(k)))%计算平衡点x%|f'(x)|<1则稳定%Exercise18%先写M文件functionf=ex4_18(a,x0,n)x=zeros(1,n);y=x;x(1)=x0;y(1)=a*x(1)+1;x(2)=y(1);plot([x(1),x(1),x(2)],[0,y(1),y(1)],'r');holdon;fori=2:ny(i)=a*x(i)+1;x(i+1)=y(i);plot([x(i),x(i),x(i+1)],[y(i-1),y(i),y(i)])endholdoff;%再执行指令>>ex4_18(0.9,1,20)>>ex4_18(-0.9,1,20)>>ex4_18(1.1,1,20)>>ex4_18(-1.1,1,20)%Exercise19clear;close;x(1)=0;y(1)=0;fork=1:3000x(k+1)=1+y(k)-1.4*x(k)^2;y(k+1)=0.3*x(k);endplot(x(1000:1500),y(1000:1500),'+g');holdonplot(x(1501:2000),y(1501:2000),'.b');plot(x(2001:2500),y(2001:2500),'*y');plot(x(2501:3001),y(2501:3001),'.r');%Exercise1x=[0410121522283440];y=[013689530];trapz(x,y)%Exercise2x=[0410121522283440];y=[013689530];diff(y)./diff(x)%Exercise3xa=-1:0.1:1;ya=0:0.1:2;[x,y]=meshgrid(xa,ya);z=x.*exp(-x.^2-y.^3);[px,py]=gradient(z,xa,ya);px%Exercise4t=0:0.01:1.5;x=log(cos(t));y=cos(t)-t.*sin(t);dydx=gradient(y,x)[x_1,id]=min(abs(x-(-1)));%找最接近x=-1的点dydx(id)%Exercise5(2)fun=inline('exp(2*x).*cos(x).^3');quadl(fun,0,2*pi)或用trapzx=linspace(0,2*pi,100);y=exp(2*x).*cos(x).^3;trapz(x,y)%Exercise5(3)fun=@(x)x.*log(x.^4).*asin(1./x.^2);quadl(fun,1,3)或用trapzx=1:0.01:3;y=feval(fun,x);trapz(x,y)Exercise5(4)fun=@(x)sin(x)./x;quadl(fun,1e-10,1)%注意由于下限为0，被积函数没有意义，用很小的1e-10代替%Exercise5(5)%参考Exercise5(4)%Exercise5(6)fun=inline('sqrt(1+r.^2.*sin(th))','r','th');dblquad(fun,0,1,0,2*pi)%Exercise5(7)首先建立84页函数dblquad2clear;fun=@(x,y)1+x+y.^2;clo=@(x)-sqrt(2*x-x.^2);dup=@(x)sqrt(2*x-x.^2);dblquad2(fun,0,2,clo,dhi,100)%Exercise6t=linspace(0,2*pi,100);x=2*cos(t);y=3*sin(t);dx=gradient(x,t);dy=gradient(y,t);f=sqrt(dx.^2+dy.^2);trapz(t,f)%Exercise7xa=-1:0.1:1;ya=0:0.1:2;[x,y]=meshgrid(xa,ya);z=x.*exp(x.^2+y.^2);[zx,zy]=gradient(z,xa,ya);f=sqrt(1+zx.^2+zy.^2);s=0;fori=2:length(xa)forj=2:length(ya)s=s+(xa(i)-xa(i-1))*(ya(j)-ya(j-1))*(f(i,j)+f(i-1,j)+f(i,j-1)+f(i-1,j-1))/4;endends%Exercise8funl=inline('-(-x).^0.2.*cos(x)');funr=inline('x.^0.2.*cos(x)');quadl(funl,-1,0)+quadl(funr,0,1)%Exercise9(以I32为例)fun=@(x)abs(sin(x));h=0.1;x=0:h:32*pi;y=feval(fun,x);t1=trapz(x,y)h=pi;x=0:h:32*pi;y=feval(fun,x);t2=trapz(x,y)%步长与周期一致，结果失真q1=quad(fun,0,32*pi)q2=quadl(fun,0,32*pi)%Exercise10(2)%先在程序编辑器，写下列函数，保存为ex5_10_2ffunctiond=ex5_10_2f(fname,a,h0,e)h=h0;d=(feval(fname,a+h)-2*feval(fname,a)+feval(fname,a-h))/(h*h);d0=d+2*e;whileabs(d-d0)>ed0=d;h0=h;h=h0/2;d=(feval(fname,a+h)-2*feval(fname,a)+feval(fname,a-h))/(h*h);end%再在指令窗口执行fun=inline('x.^2*sin(x.^2-x-2)','x');d=ex5_10_2f(fun,1.4,0.1,1e-3)%Exercise11%提示：f上升时，f'>0;f下降时，f'<0;f极值，f'=0.%Exercise12在程序编辑器，写下列函数，保存为ex5_12ffunctionI=ex5_12(fname,a,b,n)x=linspace(a,b,n+1);y=feval(fname,x);I=(b-a)/n/3*(y(1)+y(n+1)+2*sum(y(3:2:n))+4*sum(y(2:2:n)));%再在指令窗口执行ex5_12(inline('1/sqrt(2*pi)*exp(-x.^2/2)'),0,1,50)%Exercise13fun=inline('5400*v./(8.276*v.^2+2000)','v');quadl(fun,15,30)%Exercise14重心不超过凳边沿。1/2,2/3,3/4,...,n/(n+1)%Exercise15利润函数fun=inline('(p-c0+k*log(M*exp(-a*p)))*M*exp(-a*p)','p');求p使fun最大%Exercise16clear;x=-3/4:0.01:3/4;y=(3/4+x)*2.*sqrt(1-16/9.*x.^2)*9.8;P=trapz(x,y)%单位：千牛%Exercise17clear;close;fplot('17-t^(2/3)-5-2*t^(2/3)',[0,20]);grid;t=fzero('17-x^(2/3)-5-2*x^(2/3)',7)t=0:0.1:8;y=17-t.^(2/3)-5-2*t.^(2/3);trapz(t,y)-20%单位：百万元%Exercise18%曲面面积计算%Excercise1(1)fun=inline('x+y','x','y');[t,y]=ode45(fun,[0123],1)%注意由于初值为y(0)=1,[0123]中0不可缺%Excercise1(3)%令y(1)=y,y(2)=y',化为方程组%y(1)'=y(2),y(2)'=0.01*y(2)^2-2*y(1)+sin(t)%运行下列指令clear;close;fun=@(t,y)[y(2);0.01*y(2)^2-2*y(1)+sin(t)];[t,y]=ode45(fun,[05],[0;1]);plot(t,y(:,1))%Excercise1(5)%令y(1)=y,y(2)=y',化为方程组%y(1)'=y(2),y(2)'=-mu*(y(1)^2-1)*y(2)-y(1)%运行下列指令,注意参数mu的处理clear;close;fun=@(t,y,mu)[y(2);-mu*(y(1)^2-1)*y(2)-y(1)];[t,y]=ode45(fun,[020],[2;0],[],1);plot(y(:,1),y(:,2));holdon;[t,y]=ode45(fun,[020],[2;0],[],2);plot(y(:,1),y(:,2),'r');holdoff;%Excercise2roots([1105413213750])%通解A1*exp(-3*t)*cos(4*t)+A2*exp(-3*t)*sin(4*t)+A3*exp(-2*t)+A4*exp(-t)+A5*t*exp(-t)%Excercise3dfun=inline('[-1000.25*y(1)+999.75*y(2)+0.5;999.75*y(1)-1000.25*y(2)+0.5]','x','y');[x,y]=ode45(dfun,[0,50],[1;-1]);length(x)[x,y]=ode15s(dfun,[0,50],[1;-1]);length(x)%所用节点很多%所用节点很少%Excercise4clear;dfun=inline('[x(2);2*x(3)+x(1)-((1-1/82.45)*(x(1)+1/82.45))/(sqrt((x(1)+1/82.45)^2+x(3)^2))^3-(1/82.45*(x(1)-1+1/82.45))/(sqrt((x(1)+1-1/82.45)^2+x(3)^2))^3;x(4);-2*x(2)+x(3)-((1-1/82.45)*x(3))/(sqrt((x(1)+1/82.45)^2+x(3)^2))^3-(1/82.45*x(3))/(sqrt((x(1)+1-1/82.45)^2+x(3)^2))^3]','t','x');[t,x]=ode45(dfun,[024],[1.2;0;0;-1.04935371]);plot(x(:,1),x(:,3));%Excercise5%方程y'=2x+y^2,y(0)=0clear;close;fun=inline('2*x+y^2','x','y');[x,y]=ode45(fun,[01.57],0);%x的上界再增加,解会"爆炸"plot(x,y)%Excercise6clear;close;fun=@(t,x,a,b)a*x+b;[t,x]=ode45(fun,[010],0.1,[],1,1);subplot(2,4,1);plot(t,x)[t,x]=ode45(fun,[010],-0.1,[],1,1);subplot(2,4,2);plot(t,x)[t,x]=ode45(fun,[010],0.1,[],1,-1);subplot(2,4,3);plot(t,x)[t,x]=ode45(fun,[010],-0.1,[],1,-1);subplot(2,4,4);plot(t,x)[t,x]=ode45(fun,[010],0.1,[],-1,1);subplot(2,4,5);plot(t,x)[t,x]=ode45(fun,[010],-0.1,[],-1,1);subplot(2,4,6);plot(t,x)[t,x]=ode45(fun,[010],0.1,[],-1,-1);subplot(2,4,7);plot(t,x)[t,x]=ode45(fun,[010],-0.1,[],-1,-1);subplot(2,4,8);plot(t,x)%Excercise7%微分方程T'=k(c-T),T(0)=20dsolve('DT=k*(c-T)','T(0)=20','t')%得c+exp(-k*t)*(-c+20)%利用T(10)=25.2,T(20)=28.32拟合(或者解非线性方程)fun=inline('c(1)+exp(-c(2)*t)*(-c(1)+20)','c','t')lsqcurvefit(fun,[301],[1020],[25.228.32])%解得户外温度c=33,比例系数k=0.05.%Excercise8%微分方程x'/x=0.5*(1-x),x(0)=0.05fun=inline('0.5*(1-x)*x','t','x');[t,x]=ode45(fun,[010],0.05);plot(t,x)id=min(find(x>0.8));t(id)%Excercise9%微分方程组V'(t)=K(t)*V(t)^a,K'(t)=-b*K(t)%答案(1)exp(20);(2)0.353;(3)30;(4)451,0.4,9.6%Chapter7%Exercise1symsphth;a=sin(ph)*cos(th)-cos(ph)*sin(th)-sin(ph-th);simple(a)%化简后差的结果为0%Exercise2symsx;s=x^4-5*x^3+5*x^2+5*x-6;factor(s)%Exercise3symsa;A=[12;2a];iA=inv(A),[v,d]=eig(A)%Exercise4symsxy;limit((3^x+9^x)^(1/x),x,inf)s1=limit(x*y/(sqrt(x*y+1)-1),x,0);s2=limit(s1,y,0)%Exercise5symsknx;s1=symsum(k^2,k,1,n);s1=simple(s1)s2=symsum(k^(-2),k,1,inf);s2=simple(s2)s3=symsum(1/(2*n+1)/(2*x+1)^(2*n+1),n,0,inf);s3=simple(s3)%Exercise6symsxyz;s=sin(x^2*y*z);s=diff(s,x,2);s=diff(s,y,1);s=subs(s,{x,y,z},{1,1,3})%Exercise7symsx;s=log(x+sqrt(1+x^2));taylor(s,8,0,x)%Exercise8(以第四章习题9为例)先用符号运算求偏导数symsxy;f=y^3/9+3*x^2*y+9*x^2+y^2+x*y+9;fx=diff(f,x),fy=diff(f,y)根据计算结果得方程组.求解方程组[sx,sy]=solve(6*x*y+18*x+y,1/3*y^2+3*x^2+2*y+x,x,y)得四个解(0,0),(-1/3,-6),(-7/6,-7/2),(5/6,-5/2).计算Hesse矩阵fh2=[diff(fx,x),diff(fx,y);diff(fy,x),diff(fy,y)]计算eig(subs(fh2,[x,y],[0,0]))得知正定，所以是极小值点.极小值用subs(f,[x,y],[0,0])求得。同理可得(-1/3,-6)为极大值点，其它两个为鞍点。%Exercise9(以第一小题为例)symsy;f=exp(2*y)/(exp(y)+2);fi=int(f,y)s=simple(diff(fi)-f)%Exercise10symsxy;f=(x-y)^3*sin(x+2*y);Ix=simple(int(f,y,-x,x))%Exercise11(3)symsx;f=x*log(x^4)*asin(1/x);Ix=int(f,x,1,3);vpa(Ix)%Exercise12%1(3)symsx;solve(5*x^23-6*x^7+8*x^6-5*x^2)%6symsab;s=solve(a-0.7*sin(a)-0.2*cos(b),b-0.7*cos(a)+0.2*sin(b));s.a,s.b%Exercise13%1(3)dsolve('D2y-0.01*Dy^2+2*y=sin(t)','y(0)=0','Dy(0)=1','t')(解不出)%1(4)dsolve('2*D2x-5*Dx+3*x=45*exp(2*t)','x(0)=2','Dx(0)=1','t')%Exercise14%6(ii)ezplot('x^2/4+y^2/9=1')%6(vi)ezmesh('2*sin(ph)*cos(th)','2*sin(ph)*sin(th)','2*cos(ph)',[0pi/202*pi])