MAAB数学实验第二版答案胡良剑

1、数学实验答案Chapter 1Page20,ex1(5)等于exp(1 ),exp(2);exp(3),exp(4) 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)Page20, ex2(1) a, b,c的值尽管都是

2、1,但数据类型分别为数值，字符，逻辑，注意a与c相等，但他们不等于b(2) double(fun)输出的分别是字符a,b,s,(,x,)的ASCII码 Page20,ex3» r=2;p=0.5;n=12;» T=log(r)/n/log(1 +0.01 *p)Page20,ex4» x=-2:0.05:2;f=x.A4-2.Ax;» fmin,min_index=min(f)最小值最小值点编址» x(minjndex)ans =求近似根-绝对值最小的点0.6500最小值点 » f1 ,x1 Jndex=min(abs(f)f1 =0

3、.0328x1 Jndex =24» x(x1 Jndex)ans =-0.8500删去绝对值最小的点以求函数绝对值次小的» x(x1 J n dex)=;f=x.A4-2 x'. 5求另一近似根-函数绝对值次小的点点 儿、» f2,x2_index=min(abs(f)f2 =0.0630 x2 Jndex =65» x(x2_index) ans =1.2500Page20,ex5 92 99m两岫砌 74 51 58 4098 80 7 14 16 73 55 57 64 41 4 81 88 20 22 54 5663 70 4785 8

4、7 19 21 3 60 62 69 71 2886 93 25 2 9 61 68 75 52 3417 24 76 83 90 42 49 26 33 6523 5 82 89 91 48 30 32 39 6679 6 13 95 97 29 31 38 45 7210 12 94 96 78 35 37 44 46 5311 18 100 77 84 36 43 50 27 59» sum(z)» sum(diag(z)» z(:,2)/sqrt(3)» z(8,:)=z(8,:)+z(3,:)Chapter 2Page 45 ex1先在编辑器窗

5、口写下列M函数，保存为 eg2_1.mfunction xbar,s=ex2_1 (x)n=length(x);xbar=sum(x)/n;s=sqrt(sum(x42)-n*xba"2)/(n-1);例如»X=81 70 65 51 76 66 90 87 61 77；»xbar,s=ex2_1 (x)Page 45 ex2 s=log(1);n=0;while s<=100 n=n+1;s=s+log(1+n);endm=nPage 40 ex3 clear;F(1)=1;F(2)=1 ;k=2;x=0;e=1e-8; a=(1+sqrt(5)/2;wh

6、ile abs(x-a)>ek=k+1;F(k)=F(k-1)+F(k-2); x=F(k)/F(k-1); enda,x,k计算至k=21可满足精度Page 45 ex4 clear;tic;s=0;for 1=1:1000000s=s+sqrt(3)/29;ends,toctic;s=0;i=1;while i<=1000000 s=s+sqrt(3)/29;i=i+1;end s,toctic;s=O;1=1:1000000；s=sqrt(3)*sum(1 ./2.Ai);s,tocPage 45 ex5t=0:24；c=15 14 14 14 14 15 16 18 20

7、22 23 25 28 .31 32 31 29 27 25 24 22 20 18 17 16；plot(t,c)Page 45 ex6(1)x=-2:0.1:2;y=x.A2.*sin(x.A2-x-2);plot(x,y)y=inline(,xA2*sin(xA2-x-2)')；fplot(y,-2 2)(2)参数方法t=linspace(0,2*pi,100);x=2*cos(t);y=3*sin(t); plot(x,y)x=-3:0.1:3;y=x;x,y=meshgrid(x,y);z=x.A2+y.A2;surf(x,y,z) (4)x=-3:0.1:3;y=-3:0.

8、1:13;x,y=meshgrid(x,y);z=x.A4+3*x.A2+y.A2-2*x-2*y-2*x.A2.*y+6; surf(x,y,z)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

9、(05pi,100); y1 =sin(x);y2=sin(x).*sin(10*x);y3=-sin(x);plot(x5y1,x,y25xJy3)page45,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(x5y)page45,ex9clear;close;x=-2:0.1:2;y=x;x,y=meshgrid(x5y);a=0.5457；b=0.7575;p=a*exp(-0.75*y42-3.75*x.A2-1.5*x).*(x+y>1);p=p

10、+b*exp(-y.A2-6*x42).*(x+y>-1 ).*(x+yv=1);p=p+a*exp(-0.75*y.A2-3.75*x.A2+1.5*x).*(x+y<=-1);mesh(x,y,p)page45, ex10lookfor lyapunovhelp lyap» A=1 2 3；4 5 6；7 8 0；C=2 -5 -22；-5 -24 -56；-22 -56 -16；» X=lyap(A,C)X =1.0000 -1.0000 -0.0000-1.0000 2.0000 1.0000-0.0000 1.0000 7.0000Chapter 3P

11、age65 Ex1» a=1,2,3;b=2,4,3;a./b,aAb5a/b,abans = 0.5000 0.5000 1.0000ans = 221 ans =0.6552 一元方程组 x2,4,3=1,2,3的近似解 ans =0 0 00 0 00.6667 1.3333 1.0000的特解矩阵方程1,2,3x11 ,x12,x13;x21 ,x22,x23;x315x325x33=2,453Page65 Ex 2(1)» A=4 1 -1；3 2 -6；1 -5 3；b=9；-2；1；» rank(A), rank(A,b) A,b为增广矩阵 ans

12、=3ans =3可见方程组唯一解» x=Abx =2.38301.48942.0213» A=4 -3 3;3 2 -6；1 -5 3；b=-1 ；-2；1；» rank(A), rank(A,b) ans =3ans =3可见方程组唯一解» x=Ab x =-0.4706-0.2941 0» A=4 1；3 2；1 -5；b=1；1；1；» rank(A), rank(A,b) ans =2 ans =3可见方程组无解» x=Ab x =0.3311-0.1219最小二乘近似解(4)» a=251,-151 ;1

13、,2,1 ,-1 ;1,152,1;b=1 2 31;%注意 b 的写法» rank(a),rank(a,b) ans =3ans =3 rank(a)=rank(a,b)<4 说明有无穷多解» abans =1010 一个特解Page65 Ex3» a=251,-151 ；1,2,1 ,-1 ；1,152,1；b=123p；» x=null(a),xO=abx =-0.62550.6255-0.20850.4170X0 =10通解kx+xOPage65 Ex 4» x0=0.2 0.8,；a=0.99 0.05;0.01 0.95；&#

14、187; x1=a*x, x2=aA2*x, x10=aA10*x» x=xO;for i=1:1000,x=a*x;end,xx =0.83330.1667» XO=O.8 0.2'；» x=xO;for i=1:10003x=a*x;end,xx =0.83330.1667»v,e=eig(a)v =0.9806 -0.70710.1961 0.7071e =1.0000 0 0 0.9400» v(:,1)./xans =1.17671.1767成比例，说明X是最大特征值对应的特征向量Page65 Ex5用到公式(3.11)(3.

15、12)» B=6,2,1 ；2.25,1 ,0.2；3,0.2,1.8；x=25 5 20f；» C=B/diag(x)C =0.2400 0.4000 0.05000.0900 0.2000 0.01000.1200 0.0400 0.0900» A=eye(3,3)-CA =0.7600 -0.4000 -0.0500-0.0900 0.8000 -0.0100-0.1200 -0.0400 0.9100» D=17 17 17,;x=ADx =37.569625.7862 24.7690Page65 Ex 6(1)» a=4 1 -1 ;

16、3 2 -6;1 -5 3;det(a),inv(a),v,d=eig(a) ans =-94ans =0.2553 -0.0213 0.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.4158 0.6170d =-3.0527 0 00 3.6760 00 0 8.3766» a=1 1 -1 ;0 2 -1 ;-1 2 O;det(a),inv(a),v,d=eig(a) ans =2.0000 -2.0000 1.

17、0000 1.0000 -1.0000 1.00002.0000 -3.0000 2.0000-0.5773 0.5774 + O.OOOOi 0.5774 - O.OOOOi -0.5773 0.5774 0.5774-0.5774 0.5773 - O.OOOOi 0.5773 + O.OOOOi d =1.0000 0 00 1.0000 +O.OOOOi 00 0 1.0000-O.OOOOi(3)» A=5 7 6 5；7 10 8 7；6 8 10 9；579 10 A =57657 10876 8 10 9579 10» det(A),inv(A), v,d=

18、eig(A) ans =1ans =68.0000 -41.0000 -17.0000 10.0000 -41.0000 25.0000 10.0000 -6.0000-17.0000 10.0000 5.0000 -3.000010.0000 -6.0000 -3.0000 2.0000V =0.8304 0.0933 0.3963 0.3803-0.5016-0.3017 0.6149 0.5286-0.2086 0.7603 -0.2716 0.55200.1237 -0.5676 -0.6254 0.5209 d =0.0102 0 0 00 0.8431 0 00 0 3.8581

19、00 0 0 30.2887(以n=5为例)方法一(三个for)n=5;for i=1 :n, a(i,i)=5;endfor i=1 :(n-1)5a(i,i+1)=6;endfor i=1:(n-1)5a(i+1 ,i)=1 ;enda方法二(一个 for) n=5;a=zeros(n,n);a(1,1:2)=5 6;for i=2:(n-1 ),a(i,i-1 J,i+1 )=1 5 6;end a(n,n-1 n)=1 5;方法三(不用 for) n=5;a=diag(5*ones(n,1);b=diag(6*ones(n-1,1);c=diag(ones(n-1,1);a=a+zer

20、os(n-1,1 ),b;zeros(15n)+zeros(1 ,n);c,zeros(n-1,1)下歹ij计算» det(a)ans =665» inv(a)ans =0.3173 -0.5865 1.0286 -1.6241 1.9489-0.0977 0.4887 -0.8571 1.3534 -1.62410.0286 -0.1429 0.5429 -0.8571 1.0286-0.0075 0.0376 -0.1429 0.4887 -0.58650.0015-0.0075 0.0286 -0.0977 0.3173» v,d=eig(a)-0.7843

21、 -0.7843 -0.9237 0.9860 -0.92370.5546 -0.5546 -0.3771 -0.0000 0.3771-0.2614-0.2614 0.0000 -0.1643 0.00000.0924 -0.0924 0.0628 -0.0000 -0.0628-0.0218-0.0218 0.0257 0.0274 0.0257d =0.7574 0 0 0 00 9.2426 0 0 00 0 7.4495 0 00 0 0 5.0000 00 0 0 0 2.5505Page65 Ex 7» a=4 1 -1 ;3 2 -6;1 -5 3;v,d=eig(a

22、)v =0.0185 -0.9009 -0.3066-0.7693 -0.1240 -0.7248-0.6386 -0.4158 0.6170d =» det(v) ans =-0.9255 %v行列式正常，特征向量线性相关，可对角化» inv(v)*a*v 验算ans =-3.0527 0.0000 -0.00000.0000 3.6760 -0.0000-3.0527 0 00 3.6760 0-0.0000 -0.0000 8.3766»v2,d2=jordan(a)也可用 jordanv2 =0.0798 0.0076 0.91270.1886 -0.31

23、41 0.1256-0.1605 -0.2607 0.4213 特征向量不同d2 =8.3766 0 00 -3.0527 - O.OOOOi 00 0 3.6760 + O.OOOOi» v2a*v2ans =8.3766 0 0.00000.0000 -3.0527 0.00000.0000 0.0000 3.6760» v(:,1)./v2(:,2)对应相同特征值的特征向量成比例ans =2.44912.44912.4491» a=1 1 -1 ;0 2 -1 ;-1 2 O;v,d=eig(a)v =-0.5773 0.5774 + O.OOOOi 0.5

24、774 - O.OOOOi-0.5773 0.5774 0.5774-0.5774 0.5773 - O.OOOOi 0.5773 + O.OOOOid =1.0000 0 00 1.0000 + O.OOOOi 00 0 1.0000 - O.OOOOi» det(v)ans =5.0566e-028-5.1918e-017i v的行列式接近0,特征向量线性相关，不可对 角化» v,d=jordan(a) 1 -1 0V =1 1 00 1 10 0 1 jordan标准形不是对角的，所以不可对角化» A=5 7 6 5；7 10 8 7；6 8 10 9；5

25、7 9 10A =57657 10 8 76 8 10 9579 10» v,d=eig(A) v =0.8304 0.0933 0.3963 0.3803-0.5016-0.3017 0.6149 0.5286-0.2086 0.7603 -0.2716 0.55200.1237 -0.5676 -0.6254 0.5209 d =0.0102 0 0 00 0.8431 0 00 0 3.8581 0 0 0 0 30.2887 » inv(v)*A*vans =0.0102 0.0000 -0.0000 0.00000.0000 0.8431 -0.0000 -0.0

26、000-0.0000 0.0000 3.8581 -0.0000-0.0000 -0.0000 0 30.2887本题用jordan不行，原因未知(4)参考6(4)和7(1)Page65 Exercise 8只有(3)对称，且特征值全部大于零，所以是正定矩阵Page65 Exercise 9(1)» a=4 -3 1 3;2 -1 3 5；1 -1 -1 -1 ;3 -2 3 4；7 -6 -7 0» rank(a)ans =3» rank(a(1:3,:)ans =行为最大无关组» rank(a(1 2 4,:) 1,2,4 ans = 3»

27、 b=a(1 2 4,:),;c=a(3 5,:),;» bc线性表示的系数ans =0.5000 5.0000-0.5000 1.00000 -5.0000Page65 Exercise 10» a=1 -2 2；-2 -2 4；2 4 -2» v,d=eig(a)v =0.3333 0.9339 -0.12930.6667 -0.3304 -0.6681-0.6667 0.1365 -0.7327d =-7.0000 0 00 2.0000 00 0 2.0000» v'*v ans =1.0000 0.0000 0.0000 0.0000

28、1.0000 00.0000 0 1.0000 V确实是正交矩阵Page65 Exercise 11设经过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=1 0020000 0；0 0 1 0 3 0 0 0 0；1 0 -1 0 0 -3 0 0 0； 1 -1 0 0 0 0-40 0；0 -1 1 0 0 0 0 -5 0；0 1 0 0 0 0 0 0 -3； 0 0 0 1 0 -1 -1 0 0

29、；0 0 0 0 -1 -1 0 1 0；0 0 0 0 0 0 -1 -1 1；»b=20 50 0 0 0 0 0 0,; Abans =13.34536.44018.54203.3274-1.18071.6011 1.7263 0.42042.1467Page65 Exercise 12» A=1 2 3;4 5 6；7 8 0；» left=sum(eig(A), right=sum(trace(A)left =6.0000right =6应删去» left 二 prod(eig(A), right 二 det(A) 原题有错，(-1)5left

30、 =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.0853 0.1421 0.02840.1421 0.1421 0-0.0568 -0.1137 0.1705» norm(fA) f(A)范数接近 0ans =2.9536e-013Chapter 4Page84 Exercise 1roots(1 1 1)roots(3 0-4 0 2-1)(3)p=zeros(1,24);p(1 17 18 22)=5 -6 8-5；roots(p)p1=2 3

31、；p2=conv(p1, p1);p3=conv(p1, p2);p3(end)=p3(end)-4; % 原 p3 最后一个分量-4roots(p3)Page84 Exercise 2fun 二 in li neCx*log(sqrt(xA2-1 )+x)-sqrt(xA2-1 )-0.5*x');fzero(fun,2)Page84 Exercise 3fun 二 inlin eCxA4-2Ax');fplot(fun,-2 2);grid on;fzero(fun,-1 ),fzero(fun,1 ),fminbnd(fun,0.5,1.5)Page84 Exercise

32、 4 fun=inline('x*sin(1/x),x,);fplot(fun, -0.1 0.1);x=zeros(1,10);for i=1:10, x(i)=fzero(fun,(i-0,5)*0.01);end; x=x,-xPage84 Exercise 5fun 二 inlin eC9*x(1)A2+36*x A2+4*x(3)A2-36;x(1 )A2-2*x A2-20*x(3);16 *x(1 )-x(1F3-2*x(2)A2-16*x (3)/1 2:乂)；a,b,c=fsolve(fun,0 0 0)Page84 Exercise 6fun=(x)x(1 )-0.

33、7*sin(x(1 )-0.2*cos(x(2),x(2)-0.7*cos(x(1 )+0.2*sin(x);a,b,c=fsolve(fun,0.5 0.5)Page84 Exercise 7clear; 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,y15x2,y2); grid on;作图发现4个解的大致位置，然后分别求解y1 =fsolve('(x(1 )-2)A2+(x(2)-3+2*x(1 )A2-

34、5,2*(x(1 )-3)A2+(x(2)/3)A2-4',1.5,2) y2=fsolve('(x(1 )-2)A2+(x(2)-3+2*x(1 )A2-5,2*(x(1 )-3)A2+(x(2)/3)A2-4'51.8,-2) y3=fsolve('(x(1 )-2)A2+(x(2)-3+2*x(1 )A2-5,2*(x(1 )-3)A2+(x(2)/3)A2-4f ,3.5,-5) y4 二 fsolve(,(x(1 )-2F2+(x(2)-3+2*x(1 )F2-5,2*(x(1 )-3)A2+(x(2)/3)A2-4,M-4)Page84 Exerci

35、se 8clear;fun 二 i nlin e('x.A2.*si n(x.A2-x-2)');fplot(fun,-2 2);grid on;作图观察x(1)=-2;x(3)=fminbnd(fun,-1 ,-0.5);x(5)=fminbnd(funJ1,2);fun2=inline(,-x.A2.*sin(x.A2-x-2)t);x(2)=fminbnd(fun2,-2,-1);x(4)=fminbnd(fun2,-0.5,0.5);x(6)=2feval(fun,x)答案:以上x(3)是局部极小，x是局部极大，从最后一句知道x全局最小，x(2)最大。(2)clear;

36、fun=inline(,3*x.A5-20*x.A3+10,);fplot(fun,-3 3);grid on;作图观察x(1)=-3;x(3)=fminsearch(fun,2.5);fun 2=i nlin e('-(3*x.A5-20*x.A3+10)1);x(2)=fminsearch(fun2,-2.5);x(4)=3;feval(fun,x)(3)fun=inline('abs(xA3-xA2-x-2),);fplot(fun,0 3);grid on;作图观察fminbnd(fun,1.5,2.5)fun2=inline('-abs(xA3-xA2-x-2

37、),);fminbnd(fun2,0.5,1.5)Page84 Exercise 9close;x=-2:0.1:1 ;y=-7:0.1:1;x3y=meshgrid(x,y);z=y.A3/9+3*x.A2.*y+9*x.A2+y.A2+x.*y+9;mesh(x,y,z);grid on;作图观察fun=inline('x(2)A3/9+3*x(1 )A2*x(2)+9*x(1 )A2+x(2)A2+x(1 )*x(2)+9,);x=fminsearch(fun,O 0)求极小值fun2=inline(,-(x(2)A3/9+3*x(1 )A2*x(2)+9*x(1 )A2+x(2

38、)A2+x(1 )*x(2)+9)f);x=fminsearch(fun230 -5)求极大值Page84 Exercise 10clear;t=0:24;c=15 14 14 14 14 15 16 18 20 22 23 25 28 .31 32 31 29 27 25 24 22 20 18 17 16；p2=polyfit(t,c,2)p3=polyfit(t,c,3)fun 二 inlin e(*a(1 )*exp(a(2)*(t-14).A2) VaVf);a=lsqcurvefit(fun,O 0,t,c)初值可以试探f=feval(fun, a,t)norm(f-c)拟合效果p

39、lot(t,c,t,f)作图检验fun2=inline('b(1)*sin(pi/12*t4-b(2)+20,；b,'t,);原题修改 f(x)+20b=lsqcurvefit(fun2,0 O,t,c)figuref2=feval(fun2, b,t)norm(f2-c)拟合效果plot(t,c,t5f2)作图检验Page84 Exercise 11fun=inline('(1-x)*sqrt(10.52+x)-3.06*x*sqrt(1+x)*sqrt(5)');x=fzero(fun, 0, 1)Page84 Exercise 12住了颖研蟒00;N鹿屋总

40、价格y=x*0.3首付款额xO=x-y贷款总额a=(1 +r)AN*r*xO/(1 +r)AN-1)月付还款额r1 =4.05/12/100;x1=10*10000;公积金贷款a1=(1 +r1)AN*r1 *x1/(1 +r1)AN-1)x2=x0-x1商业贷款a2=(1 +r)AN*r*x2/(1 +r)AN-1)a=a1+a2Page84 Exercise 13列方程 th*RA2+(pi-2*th)*rA2-R*r*sin(th)=pi*rA2/2化简得 sin(2*th)-2*th*cos(2*th)=pi/2以下Matlab计算clear;fun= inlinefsin(2*th)

41、-2*th*cos(2*th)-pi/27th,)th=fsolve(fun,pi/4)R=20*cos(th)Page84 Exercise 14 function x=secant(fname,x0,x1 ,e)先在Editor口写M函数保存while abs(x0-x1)>e, x=x1 -(x1 - x0)*feval(fname,x1 )/(feval(fname,x1 )-feval(fname,xO);y=x*o,3首付款额xO=x1 ;x1=x;end再在指令窗口fun 二 in li ne(*x*log(sqrt(xA2-1 )+x)-sqrt(xA2-1 )-0.5*

42、/);secant(fun,1,2,1e-8)Page84 Exercise 15作系数为a,初值为x。,从第m步到第n步迭代过程的M函数:function f=ex4_15fun(a,x0,m,n)x(1)=xo； y(1)=a*x(1)+1 ；x(2)=y(1);if m<2, plot(x(1 ),x(1 )5x(2),0,y(1 ),y(1 );hold on; endfor i=2:ny(i)=a*x(i)+1; x(i+1)=y(i);if i>m, plot(x(i),x(i),x(i+1 ),y(i-1 ),y(i),y(i)； endendhold off;M脚本

43、文件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);Page84 Exercise 16设夹角L问题转化为min f=5/sin(t)+10/cos(t)取初始值pi/4,计算如下fun=(t)5/sin(t)+10/cos(t);t,f=fminsearch(fun, pi/4)0.6709 f = 20.8097Page84 Exercise

44、 17提示：x(k+2)=f(x(k)二 aA2*x(k)*(1 -x(k)*(1 -a*x(k)*(1 -x(k)计算平衡点x|f(x)|<1则稳定Page84 Exercise 18先写M文件function f=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'); hold on;for i=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)

45、,y(i)end hold off;再执行指令» ex4_18(0.9,1,20)» ex4_18(-0.9,1,20)» ex4_18(1.1,1,20)» ex4_18(-1.1,1,20)Page84 Exercise 19 clear; close; x(1)=0; y(1)=0;for k=1:3000x(k+1)=1 +y(k)-1.4*x(kF2; y(k+1)=0.3*x(k);endplot(x(1000:1500),y(1000:1500)；+g,);hold onplot(x(1501:2000),y(1501:2000)；.b&#

46、39;);plot(x(2001:2500),y(2001:2500)；*y');plot(x(2501 ：3001 ),y(2501 ：3001 )；.r)；Chapter 5Pagel 01 Exercise 1X=0 4 10 12 15 22 28 34 40；y=0 1 3 68 9 53 0；trapz(x,y)Pagel 01 Exercise 2X=0 4 10 12 15 22 28 34 40；y=0 1 3 6 8 9 5 3 0；diff(y)./diff(x)Pagel 01 Exercise 3 xa=-1:0.1:1 ;ya=0:0.1:2;x,y=mes

47、hgrid(xa,ya);z=x.*exp(-x.A2 -y.A3);px,py = gradient(z,xa,ya);pxPagel 01 Exercise 4 t=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)Pagel 01 Exercise 5Fun=inline( ' 1 /(sqrt(2*pi).*exp(-x. A2./2)1);Quadl(fun,0,1)(2) fun=inline(,exp(2*x).*

48、cos(x).A3');quadl(fun,0,2*pi)或用 trapz x=linspace(0,2*pi,100);y=exp(2*x).*cos(x)43;trapz(x,y)(3)fun 二 (x)x.*log(x44).*asi n( 17x42);quadl(fun,1,3)或用trapzx=1:0.01:3;y=feval(fun,x);trapz(x,y)(4)fun=(x)sin(x)./x;quadl(fun,1e-10,1) %注意由于下限为0,被积函数没有意义，用很小的1e-10代替(5)%参考 Exercise 5(4)(6)fun 二 in li ne(*

49、sqrt(1+r42.*si n(th)',r,；th,);dblquad(fun,0,1,0,2*pi)首先建立84页函数dblquad2 clear;fun=(x,y)1 +x+y42;clo 二(x)-sqrt(2*x-x42);dup=(x)sqrt(2*x-x.A2); dblquad2(fun,0,2,clo,dhi,100)Pagel 01 Exercise 6t=linspace(0,2*pi,100); x=2*cos(t);y=3*sin(t);dx=gradient(x,t);dy=gradient(y,t); f=sqrt(dx.A2+dy.A2);trapz(

50、t,f)Pagel 01 Exercise 7 xa=-1:0.1:1 ;ya=0:0.1:2;x,y=meshgrid(xa,ya);z=x.*exp(x.A2+y.A2);zx,zy=gradient(z,xa,ya);f=sqrt(1 +zx.A2+zy.A2);s=0;for i=2:length(xa)for j=2:length(ya)s=s+(xa(i)-xa(i-1 )*(ya(j)-ya(j-1 )*(f(i,j)+f(i-1 J)+f(iJ-1)+f(i-1 J-1)/4;end endsPagel 01 Exercise 8fun l=i nlin e('-(-x

51、).A0.2.*cos(x),);fun r=i nlin e(,x.A0.2.*cos(x)')； quadl(funl5-1,0)+quadl(funr,0,1)Pagel 01 Exercise 9 (以 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)Pagel 01 Exerci

52、se 10 先在程序编辑器，写下列函数，保存为ex5_10_2ffunction d=ex5_10_2f(fname,a,h0se) h=hO;d=(feval(fname,a+h)-2*feval(fname,a)+feval(fname,a-h)/(h*h);d0=d+2*e;while abs(d-d0)>e d0=d;h0=h;h=h0/2;d=(feval(fname,a+h)-2*feval(fname,a)+feval(fname,a-h)/(h*h); end再在指令窗口执行fun 二 i nli ne('x.A2*si n(x.A2-x-2) Vx1); d=e

53、x5_10_2f(fun,1.4,0.1,1 e-3)Pagel 01 Exercise 11 提示：f 上升时，f>O;f 下降时，f'vO;f 极值，f=0.Pagel01 Exercise 12在程序编辑器，写下列函数，保存为ex5_12ffunction I=ex5_12(fname,a5b,n) x=linspace(a,b,n+1);y=feval(fname,x); l=(b-a)/n/3*(y(1 )+y(n+1 )+2*sum(y(3:2:n)+4*sum(y(2:2:n);再在指令窗口执行ex5_12(i nlin e(,1/sqrt(2*pi)*exp(-x

54、.A2/2)'),0,1,50)Pagel01 Exercise 13 fun=inline(,5400*v./(8.276*v.A2+2000),V)； quadl(fun,15,30)Pagel 01 Exercise 14重心不超过凳边沿。1/2, 2/3, 3/4,n/(n+1)Pagel 01 Exercise15利润函数 fun=inline(,(p-cO+k*log(M*exp(-a*p)*M*exp(-a*p),'p,);求p使fun最大Pagel 01 Exercise 16 clear; x=-3/4:0.01:3/4;y=(3/4+x)*2.*sqrt(1

55、 -16/9.*x.A2)*9.8;P=trapz(x,y) % 单位:千牛Pagel 01 Exercise 17clear; close;fplot(,17-tA(2/3)-5-2*tA(2/3); 0,20); grid; t=fzero(,17-xA(2/3)-5-2*xA(2/3) 7) t=0:0.1:8;y=17-t.A(2/3)-5-2*t.A(2/3);trapz(t,y)-20 %单位:百万元Pagel 01 Exercise 18曲面面积计算Chapter 6Pagel 21 Exercise 1(1)fun=inline('x+y',x7y,);t,y=ode45(fun,0 1 2 3,1)%注意由于初值为 y(0)=1,0 1 2 3 中 0 不可缺(3)令y6=y/2)=y',化为方程组y(1)=y(2),y(2)=0.01 *y(2)A2-2*y(1)+sin(t)运行下列指令clear;close;fun=(t,y)y(2)