MATLAB课后习题答案

%Page20,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>=l,ture

(11)

答案表明：编址第2元素满足不等式(30>=20)和编址第4元素满足不等式(40>=10)

(12)

答案表明：编址第2行第1列元素满足不等式(30>=20)和编址第2行第2列元素满足不等

式(40>=10)

%Page20,ex2(l)a,b,c的值尽管都是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(l+0.01*p)T=

11.5813

%Page20,ex4»x=-2:0.05:2;f=x.A4-2.Ax;»ffmin,min_index]=min(f)fmin=

-1.3907%最小值min_index=

54%最小值点编址》x(min_index)ans=

0.6500%最小值点》[fl,xlJndex]=min(abs(D)%求近似根-绝对值最小的点fl=

0.0328xl_index=

24»x(xl_index)ans=

-0.8500»x(xl_index)=[];f=x.A4-2.Ax;%删去绝对值最小的点以求函数绝对值次小的点》

[f,2,x2_index]=min(abs(f))%求另一近似根--函数绝对值次小的点f2=

0.0630x2_index=65

»x(x2_index)ans=

1.2500

%Page20,ex5»z=magic(10)z=

929918156774515840

9880714167355576441

4818820225456637047

8587192136062697128

869325296168755234

17247683904249263365

2358289914830323966

7961395972931384572

10129496783537444653

111810077843643502759»sum(z)ans=

505505505505505505505505505505»sum(diag(z))ans

505»z(:,2)/sqrt(3)ans=

57.1577

46.1880

46.7654

50.2295

53.6936

13.8564

2.8868

3.4641

6.9282

10.3923»z(8,:)=z(8,:)+z(3,:)z=

929918156774515840

9880714167355576441

4818820225456637047

8587192136062697128

869325296168755234

17247683904249263365

2358289914830323966

83871011151198387101115119

10129496783537444653

111810077843643502759

%Page40ex1

先在编辑器窗口写下列M函数，保存为eg2」.m

functionlxbar,s]=ex2_l(x)

n=length(x);

xbar=sum(x)/n;

s=sqrt((sum(x.A2)-n*xbarA2)/(n-1));

例如

»x=[81706551766690876177];

»[xbar,sj=ex2_l(x)

xbar=

72.4000s=12.1124

%Page40ex2s=log(l);n=0;whiles<=100

n=n+1;

s=s+log(l+n);endm=n计算结果m=37

%Page40ex3

clear;

F(1)=1;F(2)=1;k=2;x=0;

e=le-8;a=(l+sqrt(5))/2;

whileabs(x-a)>e

k=k+1;F(k)=F(k-1)+F(k-2);x=F(k)/F(k-1);enda,x,k计算至k=21可满足精度

%Page40ex4clear;tic;s=0;fori=1:1000000

s=s+sqrt(3)/2Ai;ends,toctic;s=0;i=l;whilei<=1000000

s=s+sqrt(3)/2Ai;i=i+1;

end

s,toc

tic;s=0;

i=l:1000000;

s=sqrt(3)*sum(1./2.Ai);

s,toc

%Page40ex5

t=0:24;

c=[15141414141516182022232528...

313231292725242220181716];plot(t,c)

%Page40ex6

%(D

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,[-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.A2+y.A2;

surf(x,y,z)

%(4)

x=-3:0.1:3;y=-3:0.1:13;

lx,yj=meshgrid(x,y);

z=x.A4+3*x.A2+y.A2-2*x-2*y-2*x.A2.*y4-6;

surf(x,y,z)

%(5)

t=0:0.01:2*pi;

x=sin⑴;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(O,piJ00);

yl=sin(x);y2=sin(x).*sin(10*x);y3=-sin(x);

plot(x,yl,x,y2,x,y3)

%page41,ex7

x=-1.5:0.05:1.5;

y=l.l*(x>l.l)+x.*(x<=l.l).*(x>=-l.l)-l.l*(x<-l.l);

plot(x,y)

%page41,ex8

分别使用whichtrapz,typetrapz,dirC:\MATLAB7\toolbox\matlab\datafun\

%page41,ex9

clear;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.A2-3.75*x.A2-1.5*x).*(x+y>1);

p=p+b*exp(-y.八2-6*x£2).*(x+y>-l).*(x+yv=l);

p=p+a*exp(-0.75*y.A2-3.75*x.A2+1.5*x).*(x+y<=-l);

mesh(x,y,p)

%page41,exlO

lookforlyapunov

helplyap

»A=[l23;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.0000

%Chapter3

%Exercise1»a=[l,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]=[l,2,3]的近似解

ans=000000

0.66671.33331.0000%矩阵方程[l,2,3][xlI,xl2,xl3;x21,x22,x23;x31,x32,x33]=[2,4,引的特解

%Exercise2(1)

»A=[41-1;32-6;1-53];b=[9;-2;l];

»rank(A),rank([A,b])%[A,b]为增广矩阵

ans=

3ans=

3%可见方程组唯一解》x=A\bx=

2.3830

1.4894

2.0213

%Exercise2(2)

»A=[4-33;32-6;l-53];b=[-l;-2;l];

»rank(A),rank([A,b])

ans=

3ans=

3%可见方程组唯•解》x=A\bx=

-0.4706

-0.2941

0

%Exercise2(3)

»A=[41;32;1-5];b=[l;l;l];

»rank(A),rank([A,b])

ans=

2ans=

3%可见方程组无解>>x=A\bx=

0.3311

-0.1219%最小二乘近似解

%Exercise2(4)

»a=[2,l,-l,l;l,2,l,-l;l,l,2,l];b=[l23]';%注意b的写法

»rank(a),rank([a,b])

ans=

3ans=

3%rank(a)==rank([a,b])<4说明有无穷多解》a\bans=

0

0%一个特解

%Exercise3

，

»a=[2,l,-l,l;l,2,lrl;l,l,2,l];b=[l,2,31；

»x=null(a),xO=a\b

x=

-0.6255

0.6255

-0.2085

0.4170

xO=1010

%通解kx+xO

%Exercise4

»x0=[0.20,8]*;a=[0.990.05;0.010.95];

»xl=a*x,x2=aA2*x,xl0=aA10*x

»x=xO;fori=l:1000,x=a*x;end,x

x=

0.8333

0.1667»x0=[0.80.2]*;»x=xO;fori=l:1000,x=a*x;end,xx=

0.8333

0.1667»[v,e]=eig(a)v=

0.9806-0.70710.19610.7071e=1.00000

00.9400»v(:,l)./xans=

1.1767

1.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.0500

0.09000.20000.0100

0.12000.04000.0900

»A=eye(3,3)-C

A=0.7600-0.4000-0.0500-0.09000.8000-0.0100-0.1200-0.04000.9100

»D=[171717],;x=A\D

x=37.569625.786224.7690

%Exercise6(1)

»a=[41-1;32-6;1-53];det(a),inv(a),[v,d]=eig(a)ans=-94

ans=0.2553-0.02130.04260.1596-0.1383-0.22340.1809-0.2234-0.0532

v=0.0185-0.9009-0.3066-0.7693-0.1240-0.7248-0.6386-0.41580.6170

d=

-3.052700

03.67600

008.3766

%Exercise6(2)

»a=[11-l;02-1;-12O];det(a),inv(a),[v,d]=eig(a)

ans=

1

ans=2.0000-2.00001.00001.0000-1.00001.00002.0000-3.00002.0000

v=-0.57730.5774+O.OOOOi0.5774-O.OOOOi-0.57730.57740.5774-0.57740.5773-O.OOOOi

0.5773+O.OOOOi

d=

1.000000

01.0000+O.OOOOi0

001.0000-O.OOOOi

%Exercise6(3)

»A=[5765;71087;68109;57910]

A=

5765

71087

68109

57910

»det(A),inv(A),[v,d]=eig(A)ans=1

ans=68.0000-41.0000-17.000010.0000-41.000025.000010.0000-6.0000-17.000010.0000

5.0000-3.000010.0000-6.0000-3.00002.0000

v=0.83040.09330.39630.3803-0.5016-0.30170.61490.5286-0.20860.7603-0.27160.5520

0.1237-0.5676-0.62540.5209

d=0.010200000.843100003.8581000030.2887

%Exercise6(4)(以n=5为例)

%关键是矩阵的定义

%方法一(三个for)

n=5;

fori=l:n,a(i,i)=5;end

fori=1:(n-1),a(iJ+1)=6;end

fori=I:(n-l),a(i+l,i)=l;end

a

%方法二(一个for)

n=5;a=zeros(n,n);

a(l,l:2)=[56];

fori=2:(n-1),a(i,[i-1,i,i+1])=[156];end

a(n,[n-ln])=[l5];

a

%方法三(不用for)

n=5;a=diag(5*ones(n,l));

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-l,l)]

%下列计算

»det(a)

ans=665

»inv(a)

ans=0.3173-0.58651.0286-1.62411.9489-0.09770.4887-0.85711.3534-1.6241

0.0286-0.14290.5429-0.85711.0286

-0.00750.0376-0.14290.4887-0.5865

0.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.0257

0.02740.0257

d=

0.75740000

09.2426000

007.449500

0005.00000

00002.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.6170

d=

-3.052700

03.67600

008.3766

»det(v)ans=

-0.9255%v行列式正常，特征向量线性相关，可对角化》出丫")*@*丫％验算ans=

-3.05270.0000-0.0000

0.00003.6760-0.0000

-0.0000-0.00008.3766

»[v2,d2]=jordan(a)%也可用jordan

v2=0.07980.00760.91270.1886-0.31410.1256

-0.1605-0.26070.4213%特征向量不同d2=8.3766000-3.0527-O.OOOOi0

003.6760+O.OOOOi»v2\a*v2ans=

8.376600.00000.0000-3.05270.00000.00000.00003.6760

»v(:,l)./v2(:,2)%对应相同特征值的特征向量成比例

ans=2.44912.44912.4491

%Exercise7(2)

»a=[l120];[v,d]=eig(a)

v=

-0.57730.5774+O.OOOOi0.5774-O.OOOOi

-0.57730.57740.5774

-0.57740.5773-O.OOOOi0.5773+O.OOOOi

d=

1.00000001.0000+O.OOOOi0001.0000-O.OOOOi

»det(v)

ans=

-5.0566e-028-5.1918e-017i%v的行歹ij式接近0,特征向量线性相关，不可对角化

»[v,d]=jordan(a)

101

100

1-10

d=

110

011

001%jordan标准形不是对角的，所以不可对角化

%Exercise7(3)

»A=[5765;71087;68109;57910]

A=

5765

71087

68109

57910»[v,d]=eig(A)

v=0.83040.09330.39630.3803-0.5016-0.30170.61490.5286-0.20860.7603-0.27160.5520

0.1237-0.5676-0.62540.5209

d=

0.0102000

00.843100

003.85810

00030.2887

»inv(v)*A*v

ans=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(l:3,:))ans=

2»rank(a([124],:))%1,2,4行为最大无关组ans=

3»b=a([l24],:)*;c=a([35],:)';»b\c%线性表示的系数ans=

0.50005.0000

-0.50001.0000

0-5.0000

%Exercise10

»a=ll-22;-2-24;24-2]

»[v,d]=eig(a)

v=

0.33330.9339-0.1293

0.6667-0.3304-0.6681

-0.66670.1365-0.7327d=-7.00000002.00000

002.0000

»v**v

ans=

1.00000.00000.0000

0.00001.00000

0.000001.0000%v确实是正交矩阵

%Exercise11

%设经过6个电阻的电流分别为il,…,i6.列方程组如下

%20-2il=a;5-3i2=c;a-3i3=c;a-4i4=b;c-5i5=b;b-3i6=0;

%il=i3+i4;i5=i2+i3;i6=i4+i5;

%计算如下

»A=[100200000;001030000;l0-100-3000;1-10000-400;

0-110000-50;01000000-3;00010-1-100;0000-l-1010;

000000-1-11];

»b=[2050000000],;A\b

ans=

13.34536.44018.54203.3274-1.18071.60111.72630.42042.1467

%Exercise12

»A=[l23;456;780];

»left=sum(eig(A)),right=sum(trace(A))

left=

6.0000right=

»left=prod(eig(A)),right=det(A)%原题有错,(-1)八n应删去left=

27.0000right=

27»fA=(A-p(l)*eye(3,3))*(A-p(2)*eye(3,3))*(A-p(3)*eye(3,3))fA=

1.0e-012*

0.08530.14210.0284

0.14210.14210

-0.0568-0.11370.1705»norm(fA)%f(A)范数接近Oans=

2.9536e-013

%Exercise1(1)

roots([l11J)

%Exercise1(2)

roots([30-402-1])

%Exercise1(3)

p=zeros(l,24);

p([l171822])=[5-68-5];

roots(p)

%Exercise1(4)

pl=[23];

p2=conv(p1,pl);

p3=conv(pl,p2);

p3(end)=p3(end)-4;%原p3最后一个分量-4

roots(p3)

%Exercise2

fun=inline('x*log(sqrt(x八2-l)+x)-sqrt(xA2-l)-0.5*x');

fzero(fun,2)

%Exercise3

fun=inline(,xA4-2Ax,);

fplot(fun,[-22]);gridon;

fzero(fun,-1),fzero(fun,1),fminbnd(fun,0.5,1.5)

%Exercise4

fun=inline('x*sin(l/x)','x');

fplot(fun,[-0.10.1]);

x=zeros(l,10);fori=l:10,x(i)=fzero(fun,(i-0.5)*0.01);end;

x=[x,-x]

%Exercise5

fun二inline(19*x(1)八2+36*x(2)八2+4*x(3)八2-36;x(l)八2-2*x(2)八2-20*x(3);16*x(1)-x(1)A3-2*x(2)A

2-16*x(3)A2],；x,);

[a,b,c]=fsolve(fun,[000])

%Exercise6

fun=®(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])

%Exercise7

clear;close;t=0:pi/100:2*pi;

xl=2+sqrt(5)*cos(t);y1=3-2*x1+sqrt(5)*sin(t);

x2=3+sqrt(2)*cos(t);y2=6*sin(t);

plot(xl,yl,x2,y2);gridon;%作图发现4个解的大致位置，然后分别求解

yl=fsolveC[(x(1)-2)A2+(x(2)-3+2*x(1))A2-5,2*(x(1)-3)A2+(x(2)/3)A2-4],,[1.5,2])

y2=fsolve(,[(x(l)-2)A2+(x(2)-3+2*x(l))A2-5,2*(x(l)-3)A2+(x(2)/3)A2-4]',[1.8,-2])

y3=fsolveC[(x(1)-2)A2+(x(2)-3+2*x(1))A2-5,2*(x(1)-3)A2+(x(2)/3)A2-4],,[3.5,-5])

y4=fsolveC[(x(1)-2)A2+(x(2)-3+2*x(1))A2-5,2*(x(1)-3)A2+(x(2)/3)A2-4J;[4,-4J)

%Exercise8(1)

clear;

fun=inline(,x.A2.*sin(x.A2-x-2),);

fplot(fun,[-22]);gridon;%作图观察

x(l)=-2;

x(3)=fminbnd(fun,-1,-0.5);

x(5)=fminbnd(fun,1,2);

fun2=inline(,-x.A2.*sin(x.A2-x-2)*);

x(2)=fminbnd(fun2,-2,-1);

x(4)=fminbnd(fun2,-0.5,0.5);

x(6)=2

feval(fun,x)

%答案：以上x(l)(3)(5)是局部极小，x(2)(4)(6)是局部极大，从最后一句知道x(l)全局最小，

x(2)最大。

%Exercise8(2)

clear;

fun=inline(,3*x.A5-20*x.A3+10,);

fplot(fun,[-33]);gridon;%作图观察

x(l)=-3;

x(3)=fminsearch(fun,2.5);

fun2=inline('-(3*x.A5-20*x.A3+10),);

x(2)=fminsearch(fun2,-2.5);

x(4)=3;

feval(fun,x)

%Exercise8(3)

fun=inline('abs(xA3-x八2-x-2)');

fplot(fun,[03]);gridon;%作图观察

fminbnd(fun,1.5,2.5)

fun2=in1ine('-abs(x八3-x八2-x-2)');

fminbnd(fun2,0.5,1.5)

%Exercise9

close;

x=-2:0.1:l;y=-7:0.1:l;

[x,yj=meshgrid(x,y);

z=y.八3/9+3*x.八2.*y+9*x.A2+y.八2+x.*y+9;

mesh(x,y,z);gridon;%作图观察

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,[O0])%求极小值

fun2=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(fun2,[0-5])%求极大值

%Exercise10

clear;t=0:24;

c=[15141414141516182022232528...

313231292725242220181716J;p2=polyfit(t,c,2)p3=polyfit(t,c,3)

fun=inline(,a(l)*exp(a(2)*(t-14).A2),,*a,,,t,);a=lsqcurvefit(fun,[OO],t,c)%初值可以试探

fi=feval(fun,a,t)norm(f-c)%拟合效果plot(t,c,t,f)%作图检验

fun2=inline(b(l)*sin(pi/12*t+b(2))+20/b；t);%原题修改f(x)+20b=lsqcurvefit(fun2,[001,t,c)

figuref2=feval(fun2,b,t)norm(f2-c)%拟合效果plot(t,c,t,f2)%作图检验

%Exercise11fun=inline('(l-x)*sqrt(10.52+x)-3.06*x*sqrt(l+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=(l+r)八N*r*x0/((l+r)八N-1)%月付还款额r1=4.05/12/100;x1=10*10000;%

公积金贷款a1=(1+r1)AN*r1*x1/((1+r1)AN-1)x2=x0-xl%商业贷款

a2=(l+r)AN*r*x2/((l+r)AN-l)a=a1+a2

%Exercise13%列方程th*R八2+(pi-2*lh)*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,xO,x1,e)while

abs(x0-xl)>e,

x=x1-(x1-xO)*feval(fname,x1)/(feval(fname,x1)-feval(fname,xO));

x0=xl;xl=x;end%再在指令窗口fun=inline('x*log(sqrt(xA2-l)+x)-sqrt(xA2-l)-0.5*x');

secant(fun,1,2,1e-8)

%Exercise15

%作系数为a,初值为xo,从第m步到第n步迭代过程的M函数：

functionf=ex4_l5fun(a,x0,m,n)

x(l)=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;end

fori=2:n

y(i)=a*x(i)+l;x(i+l)=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,l,1.20);subplot(2,2,2);ex4_l5fun(-0.9,l,l,20);

subplot(2,2,3);ex4_15fun(l.l,l,l,20);subplot(2,2,4)；ex4_15fun(-l.l,l,l,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))=aA2*x(k)*(1-x(k))*(1-a*x(k)*(1-x(k)))

%计算平衡点x

%lf(x)l<l则稳定

%Exercise18

%先写M文件

functionf=ex4_18(a,x0,n)

x=zeros(l,n);y=x;

x(l)=x0;

y(l)=a*x(l)+l;

x(2)=y(l);

plot([x(1),x(1),x(2)],[0,y(1),y(1

holdon;

fori=2:n

y(i)=a*x(i)+l;

x(i+l)=y(i);

plot([x(i),x(i),x(i+l)],[y(i-l),y(i),y(i)])endholdoff;%再执行指令：>>ex4_18(0.9,l,20)»

ex4_l8(-0.9,l,20)»ex4_18(l.l,l,20)»ex4_18(-1.1,1,20)

%Exercise19clear;close;x(l)=0;y(l)=0;fork=1:3000

x(k+l>l+y(k)-1.4*x(k)A2;y(k+l)=0.3*x(k);endplot(x(1000:1500),y(1000:1500)；+g*);holdon

plot(x(l501:2000),y(l501:2000),'.b');plot(x(2001:2500),y(2001:2500)；*y);

plot(x(2501:3001),y(2501:3001)；.r');

%Exercise1

x=[0410121522283440];

y=[013689530];

trapz(x,y)

%Exercise2

x=[0410121522283440];

y=[013689530];

diff(y)./diff(x)

%Exercise3

xa=-1:0.1:1;ya=0:0.1:2;

[x,y]=meshgrid(xa,ya);

z=x.*exp(-x.A2-y.A3);

[px,py]=gradient(z,xa,ya);

px

%Exercise4

t=0:0.01:1.5;

x=log(cos(t));

y=cos(t)-t.*sin(t);

dydx=gradient(y,x)

[x」,id]=min(abs(x・(-l)));%找最接近x=-l的点

dydx(id)

%Exercise5(2)

fun=inline('exp(2*x).*cos(x).A3,);

quadl(fun,0,2*pi)

或用ir叩z

x=linspace(0,2*pi,100);

y=exp(2*x).*cos(x).A3;

trapz(x,y)

%Exercise5(3)

fun=®(x)x.*log(x.A4).*asin(1./x.A2);

quadl(fun,l,3)

或用trapz

x=l:0.01:3;

y=feval(fun,x)；

trapz(x,y)

%Exercise5(4)

fun=@(x)sin(x)./x;

4皿出什皿1。101)％注意由于下限为0,被积函数没有意义，用很小的le-10代替

%Exercise5(5)

%参考Exercise5(4)

%Exercise5(6)

fun=inline('sqrt(l+r.八2.*sin(th))；T;th');

dblquad(fun,0,1,0,2*pi)

%Exercise5(7)

首先建立84页函数dblquad2

clear;

fun=®(x,y)l+x+y.A2;

clo=@(x)-sqrt(2*x-x.A2);

dup=®(x)sqrt(2*x-x.A2);

dblquad2(fun,0,2,clo,dhi,100)

%Exercise6

t=linspace(0,2*pi,l00);

x=2*cos(t);y=3*sin(t);

dx=gradient(x,t);dy=gradient(y,t);

f=sqrt(dx.A2+dy.A2);

trapz(t,f)

%Exercise7

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;

fori=2:length(xa)

forj=2:length(ya)

s=s+(xa(i)-xa(i・l))*(ya(j)・yag))*(f(i,j)+f(i-l,j)+f(i,j-l)+f(i-lJ-l))/4;

endends

%Exercise8funl=inline(-(-x).A0.2.*cos(x),);

funr=inline('x.A0.2.*cos(x)');

quadl(funl,-1,0)+quadl(funr,0,1)

%Exercise9(以132为例)

fun=@(x)abs(sin(x));

h=0.1;x=0:h:32*pi;y=feval(fun,x);tl=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)

%Exercise1(X2)

%先在程序编辑器，写下列函数，保存为ex5」0_2f

functiond=ex5_10_2f(fname,a,h0,e)

h=h0;d=(feval(fname,a4-h)-2*feval(fname,a)+feval(fname,a-h))/(h:f:h);

d0=d+2*e;

whileabs(d-dO)>e

d0=d;h0=h;h=h0/2;

d=(feval(fname,a+h)-2*feval(fname,a)4-feval(fname,a-h))/(h*h);end%再在指令窗口执行

fun=inline('x.A2*sin(x.A2-x-2),,,x,);d=ex5_l0_2f(fun,1.4,0.1,1e-3)

%Exercise11

%提示：f上升时，F>O;f下降时，Fv0;f极值，f=0.

%Exercise12

在程序编辑器，写下列函数，保存为ex5_l2f

functionI=ex5_12(fname,a,b,n)

x=linspace(a,b,n+1);

y=feval(fname,x);

I=(b-a)/n/3*(y(l)+y(n+l)+2*sum(y(3:2:n))+4*sum(y(2:2:n)));

%再在指令窗口执行

ex5_12(inline(,l/sqrt(2*pi)*exp(-x.A2/2)'),0J>50)

%Exercise13

fun=inline(,5400*v./(8.276*v.A2+2000),；v');

quadl(fun,15,30)

%Exercise14

重心不超过凳边沿。1/2,2/3,3/4,…,n/(n+l)

%Exercisel5

利润函数fun=inlineC(p-cO+k*1og(M*exp(-a*p)))*M*exp(-a*p)?p);

求p使fun最大

%Exercise16

clear;x=-3/4:0.01:3/4;

y=(3/4+x)*2*sqrt(l-16/9.*x.A2)*9.8;

P=trapz(x,y)%单位：千牛

%Exercise17

clear;close;

fp!ot(,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=O:O.l:8;y=l7-t.A(2/3)-5-2*t.A(2/3);

trapz(t,y)-20%单位：百万元

%Exercise18

%曲面面积计算

%Excercise1(1)

fun=inline(,x+y,,,x,,y);

[t,y]=ode45(fun,[0123],1)%注意由于初值为y(0)=l,[0123]中0不可缺

%Excercise1(3)

%令y(D=y,y(2)=y',化为方程组

%y(l),=y(2),y(2),=0.01*y(2)A2-2*y(l)+sin(t)

%运行下列指令

clear;close;

fun=®(t,y)[y(2);0.01*y(2)A2-2*y(1)+sin(t)];

[t,y]=ode45(fun,[05],[0;l]);

plot(t,y(:,l))

%Excercise1(5)

%令y(l)=yM2)=y',化为方程组

%y(1)'=y(2),y(2)'=-mu*(y(1)A2-1)*y(2)-y(1)

%运行下列指令，注意参数mu的处理

clear;close;

fun=@(t,y,mu)[y(2);-mu*(y(1)A2-1)*y(2)-y(1)];

[t,y]=ode45(fun,[020],[2;0],[],l);

plot(y(:,l),y(:,2));holdon;

[t,y]=ode45(fun,[020],[2;0],[],2);

plot(y(:,l),y(:,2),'r');holdoff;

%Excercise2

roots([l105413213750])

%通解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)

%Excercise3

dfun=inline('[-1000.25*y(l)+999.75*y(2)+0.5;999.75*y(l)-1000.25*y(2)+0.5]'；xVy,);

[x,y]=ode45(dfun,[0,50],[1;-l]);length(x)

%所用节点很多

[x,y]=ode15s(dfun,[0,50],[l;-l]);length(x)

%所用节点很少

%Excercise4

clear;

dfun=inline(Ix(2);2*x(3)+x(l)-((l-l/82.45)*(x(1)+l/82.45))/(sqrt((x(1)+1/82.45)A2+x(3)A2))A3-(

l/82.45*(x(l)-l+l/82.45))/(sqrt((x(l)+l-l/82.45)A2+x(3)A2))A3;

x(4);-2*x(2)+x(3)-((l-l/82.45)*x(3))/(sqrt((x(l)+l/82.45)A2+x(3)A2))A3-(l/82.45*x(3))/(sqrt((x(l

)+Ll/82.45)A2+x(3)A2))A3]'；tVx,);

[t,x]=ode45(dfun,[024],[1.2;0;0;-1.04935371]);

plot(x(:,l),x(：3))；

%Excercise5

%方程y'=2x+yA2,y(0)=0

clear;close;

fun=inline(,2*x+yA2,,,x7y,);

[x,y]=ode45(fun,[01.57],0);%x的上界再增加，解会"爆炸"

plot(x,y)

%Excercise6

clear;close;

fun=@(t,x,a,b)a*x+b;

[t,x]=ode45(fun,[010],0.1,[],lJ)；

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,[],-l,D；

subplot(2,4,5);plot(t,x)

[t,xl=ode45(fun,[01

subplot(2,4,6);plot(t,x)

[t,x]=ode45(fun,[010],0.1,[],-1,-1);

subplot(2,4,7);plot(t,x)

[