MATLAB数学实验第二版答案

数学试验答案

Chapter1

Page20,exl

(5)等于[exp⑴,exp⑵;exp⑶,exp⑷]

(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元素满足不等式(3。>=20)和编址第4元素满足

不等式(40>=10)

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

行第2列元素满足不等式(40>=10)

Page20,ex2

⑴a,b,c的值尽管都是1,但数据类型分别为数值，字符，逻辑，留

意a与c相等，但他们不等于b

⑵doubleKun)输出的分别是字符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(min_index)

ans=

0.6500最小值点

»[f1,xl_index]=min(abs(f))求近似根一确定值最小的点

fl=

0.0328

xl_index=

24

>>x(xl_index)

ans=

-0.8500

»x(xl_index)=[];f=x.A4-2.Ax;删去确定值最小的点以求函数确定值

次小的点

»[f2,x2_index]=min(abs(f))求另一近似根一函数确定值次小的点

f2=

0.0630

x2_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)

>>sum(diag(z))

>>z(:,2)/sqrt(3)

>>z(8,:)=z(8,:)+z(3,:)

Chapter2

Page45exl

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

function[xbar,s]=ex2_l(x)

n=length(x);

xbar=sum(x)/n;

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

例如

»x=[81706551766690876177];

>>[xbar,s]=ex2_l(x)

Page45ex2

s=log(l);n=0;

whiles<=100

n=n+l;

s=s+log(l+n);

end

m=n

Page40ex3

clear;

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

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

whileabs(x-a)>e

k=k+l;F(k)=F(k-l)+F(k-2);x=F(k)/F(k-l);

end

a,x,k

计算至k=21可满足精度

Page45ex4

clear;tic;s=0;

fori=l:1000000

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

end

s,toc

tic;s=0;i=l;

whilei<=1000000

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

end

s,toc

tic;s=0;

i=l:1000000;

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

s,toc

Page45ex5

t=0:24;

c=[15141414141516182022232528...

313231292725242220181716];

plot(t,c)

Page45ex6

(1)

x=-2:0.l:2;y=x.A2.*sin(x.A2-x-2);plot(x,y)

y=inline('xA2*sin(xA2-x-2),);fplot(y,[-22])

⑵参数方法

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

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

⑶

x=-3:0,l:3;y=x;

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

z=x.A2+y.A2;

surf(x,y,z)

(4)

x=-3:0.l:3;y=-3:0.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(O,pi,100);

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

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

page45,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)

page45,ex9

clear;close;

x=-2:0,l: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>l);

p=p+b*exp(-y.A2-6*x.A2).*(x+y>-l).*(x+y<=l);

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

mesh(x,y,p)

page45,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

Page65Exl

>>a=[l,2,3];b=[2,4,3];a./b,a.\b,a/b,a\b

ans=

0.50000.50001.0000

ans=

221

ans=

0.6552一元方程组x[2,4,3]=[l,2,3]的近似解

ans=

000

000

0.66671.33331.0000

矩阵方程[I,2,3][xll,xl2,xl3;x21,x22,x23;x31,x32,x33]=[2,4,3]的

特解

Page65Ex2

⑴

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

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

ans=

3

ans=

3可见方程组唯一解

>>x=A\b

x=

2.3830

1.4894

2.0213

(2)

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

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

ans=

3

ans=

3可见方程组唯一解

>>x=A\b

x=

-0.4706

-0.2941

0

⑶

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

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

ans=

2

ans=

3可见方程组无解

>>x=A\b

x=

0.3311

-0.1219最小二乘近似解

(4)

»a=[2,l,-l,l;l,2,l,-l；l,l,2,l];b=[l23了;％留意b的写法

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

ans=

3

ans=

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

>>a\b

ans=

1

0

1

0一个特解

Page65Ex3

»a=[2,l,-l,l;l,2,l,-l;l,l,2,l];b=[l,2,3],;

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

x=

-0.6255

0.6255

-0.2085

0.4170

xO=

1

0

1

0

通解kx+xO

Page65Ex4

»x0=[0.20.8],;a=[0,990.05;0.010.95];

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

>>x=xO;fori=1:1000,x=a*x;end,x

x=

0.8333

0.1667

»x0=[0.80.2],;

>>x=x0;fori=1:1000,x=a*x;end,x

x=

0.8333

0.1667

»[v,e]=eig(a)

v=

0.9806-0.7071

0.19610.7071

e=

1.00000

00.9400

»v(:,l)./x

ans=

1.1767

1.1767成比例，说明x是最大特征值对应的特征向量

Page65Ex5

用到公式(3.11)(3.12)

»B=[6,2,l;2,25,l,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.5696

25.7862

24.7690

Page65Ex6

⑴

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

ans=

-94

ans=

0.2553-0.02130.0426

0.1596-0.1383-0.2234

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

(2)

>>a=[l1-l;02-1;-120];det(a),inv(a),[v,d]=eig(a)

ans=

1

ans=

2.0000-2.00001.0000

1.0000-1.00001.0000

2.0000-3.00002.0000

v=

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

-0.57730.57740.5774

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

d=

1.000000

01.0000+O.OOOOi0

001.0000-O.OOOOi

(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.00005.0000-3.0000

10.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.0102000

00.843100

003.85810

00030.2887

⑷(以n=5为例)

方法一(三个for)

n=5;

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

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

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

a

方法二(一个for)

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

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

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

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

a

方法三(不用for)

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

b=diag(6*ones(n-l,1))；

c=diag(ones(n-1,1));

a=a+[zeros(n-l,l),b;zeros(l,n)]+[zeros(l,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.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.9237

0.5546-0.5546-0.3771-0.00000.3771

-0.2614-0.26140.0000-0.16430.0000

0.0924-0.09240.0628-0.0000-0.0628

-0.0218-0.02180.02570.02740.0257

d=

0.75740000

09.2426000

007.449500

0005.00000

00002.5505

Page65Ex7

(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行列式正常，特征向量线性相关，可对角化

>>inv(v)*a*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.9127

0.1886-0.31410.1256

-0.1605-0.26070.4213特征向量不同

d2=

8.376600

0-3.0527-O.OOOOi0

003.6760+O.OOOOi

>>v2\a*v2

ans=

8.376600.0000

0.0000-3.05270.0000

0.00000.00003.6760

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

ans=

2.4491

2.4491

2.4491

⑵

>>a=[l1-l;02-1;-120];[v,d]=eig(a)

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

-0.57730.57740.5774

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

d=

1.000000

01.0000+O.OOOOi0

001.0000-O.OOOOi

>>det(v)

ans=

-5.0566e-028-5.1918e-017iv的行列式接近0,特征向量线性相

关，不行对角化

>>[v,d]=jordan(a)

v=

101

100

1-10

d=

110

011

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

⑶

»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.0000

0.00000.8431-0.0000-0.0000

-0.00000.00003.8581-0.0000

-0.0000-0.0000030.2887

本题用jordan不行，缘由未知

(4)

参考6⑷和7⑴

Page65Exercise8

只有⑶对称，且特征值全部大于零，所以是正定矩阵.

Page65Exercise9

⑴

»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([l24],:))1,2,4行为最大无关组

ans=

3

»b=a([l24],:)';c=a([35],:)1;

»b\c线性表示的系数

ans=

0.50005.0000

-0.50001.0000

0-5.0000

Page65Exercise10

»a=[l-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.7327

d=

-7.000000

02.00000

002.0000

>>v'*v

ans=

1.00000.00000.0000

0.00001.00000

0.000001.0000v的确是正交矩阵

Page65Exercise11

设经过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=[l00200000;001030000;10-100-3000;1-10

000-400;

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

-1010;

000000-1-11];

»b=[2050000000],;A\b

ans=

13.3453

6.4401

8.5420

3.3274

-1.1807

1.6011

1.7263

0.4204

2.1467

Page65Exercise12

»A=[l23;456;780];

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

left=

6.0000

right=

6

»left=prod(eig(A)),right=det(A)原题有错，(-l)An应删去

left=

27.0000

right=

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)范数接近0

ans=

2.9536e-013

Chapter4

Page84Exercise1

(1)

roots([l11])

(2)

roots([30-402-1])

⑶

p=zeros(l,24);

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

roots(p)

(4)

pl=[23];

p2=conv(pl,pl);

p3=conv(pl,p2);

p3(end)=p3(end)-4;%原p3最终一个重量-4

roots(p3)

Page84Exercise2

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

fzero(fun,2)

Page84Exercise3

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

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

fzero(fun,-l),fzero(fun,1),fminbnd(fun,0.5,1-5)

Page84Exercise4

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]

Page84Exercise5

fun=inline('[9*x(l)A2+36*x(2)A2+4*x(3)A2-36;x(l)A2-2*x(2)A2-20*x

(3);16*x(l)-x(l)A3-2*x(2)A2-16*x(3)A2]','x');

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

Page84Exercise6

fun=@(x)[x(1)-0.7*sin(x(1))-0.2*cos(x(2)),x(2)-0.7*cos(x(l))+0.2*sin(

x(2))];

[a,b,c]=fsolve(fun,[0.50,5])

Page84Exercise7

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

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

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

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

求解

yl=fsolve('[(x(l)-2)A2+(x(2)-3+2*x(l))A2-5,2*(x(l)-3)A2+(x(2)/3)A2-

4]1,[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,[1.8,-2])

y3=fsolve('[(x(l)-2)A2+(x(2)-3+2*x(l))A2-5,2*(x(l)-3)A2+(x(2)/3)A2-

4]1,[3.5,-5])

y4=fsolve('[(x(l)-2)A2+(x(2)-3+2*x(l))A2-5,2*(x(l)-3)A2+(x(2)/3)A2-

4],,[4,-4])

Page84Exercise8

(1)

clear;

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

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

x(l)=-2;

x(3)=fminbnd(fun,-l,-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⑴全局最小，x(2)最大。

(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)

(3)

fun=inline(,abs(xA3-xA2-x-2),)；

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

fminbnd(fun,1.5,2.5)

fun2=inline('-abs(xA3-xA2-x-2),);

fminbnd(fun2,0.5,1.5)

Page84Exercise9

close;

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

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

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

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

fun=inline('x(2)A3/9+3*x(l)A2*x(2)+9*x(l)A2+x(2)A2+x(l)*x(2)+9,);

x=fminsearch(fun,[00])求微小值

fun2=inline('-(x(2)A3/9+3*x(l)A2*x(2)+9*x(l)A2+x(2)A2+x(l)*x(2)+

9)1)；

x=fminsearch(fun2,[0-5])求极大值

Page84Exercise10

clear;t=0:24;

c=[15141414141516182022232528...

313231292725242220181716];

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,[00],t,c)初值可以摸索

f=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)+20

b=lsqcurvefit(fun2,[0O],t,c)

figure

f2=feval(fun2,b,t)

norm(f2-c)拟合效果

plot(t,c,t,f2)作图检验

Page84Exercise11

fun=inline('(1-x)*sqrt(l0.52+x)-3.06*x*sqrt(l+x)*sqrt(5)');

x=fzero(fun,0,1)

Page84Exercise12

r=5.04/12/100;N=20*12;

x=7500*180房屋总价格

y=x*0.3首付款额

xO=x-y贷款总额

a=(l+r)AN*r*xO/((l+r)AN-l)月付还款额

rl=4.05/12/100;xl=10*10000;公积金贷款

al=(l+rl)AN*rl*xl/((l+rl)AN-l)

x2=x0-xl商业贷款

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

a=al+a2

Page84Exercise13

列方程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=inline('sin(2*th)-2*th*cos(2*th)-pi/2','th')

th=fsolve(fun,pi/4)

R=20*cos(th)

Page84Exercise14

先在Editor窗口写M函数保存

functionx=secant(fname,xO,xl,e)

whileabs(xO-xl)>e,

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

xO=xl;xl=x;

end

再在指令窗口

fun=inline('x*log(sqrt(xA2-l)+x)-sqrt(xA2-l)-0.5*x,)；

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

Page84Exercise15

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

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

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

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)]);end

end

holdoff;

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(-l.1,1,1,20);

Page84Exercise16

设夹角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.6709

f=

20.8097

Page84Exercise17

提示：x(k+2)=f(x(k))=aA2*x(k)*(l-x(k))*(1-a*x(k)*(1-x(k)))

计算平衡点x

|f'(x)|〈l则稳定

Page84Exercise18

先写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)],'r');

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

end

holdoff;

再执行指令

»ex4_18(0.9,l,20)

»ex4_18(-0.9,l,20)

»ex4_18(l.l,l,20)

»ex4_18(-l.1,1,20)

Page84Exercise19

clear;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);

end

plot(x(1000:1500),y(1000:1500),'+g');holdon

plot(x(1501:2000),y(1501:2000),'.b');

plot(x(2023:2500),y(2023:2500),'*y')；

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

Chapter5

PagelOlExercise1

x=[0410121522283440];

y=[013689530];

trapz(x,y)

PagelOlExercise2

x=[0410121522283440];

y=[013689530];

diff(y)./diff(x)

PagelOlExercise3

xa=-l:0.l:l;ya=0:0.1:2;

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

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

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

px

PagelOlExercise4

t=0:0.01:1.5;

x=log(cos(t));

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

dydx=gradient(y,x)

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

dydx(id)

PagelOlExercise5

(1)

Fun=inline(4l/(sqrt(2*pi)).*exp(-x.A2./2),);

Quadl(fun,0,1)

(2)

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

quadl(fun,0,2*pi)

或用trapz

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

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

trapz(x,y)

(3)

fun=@(x)x.*log(x.A4).*asin(l./x.A2);

quadl(fun,l,3)

或用trapz

x=l:0.01:3;

y=feval(fun,x);

trapz(x,y)

(4)

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

口遁叫g11,16-10,1)%留意由于下限为0,被积函数没有意义，用很小的

le-10代替

(5)

%参考Exercise5(4)

(6)

fun=inline('sqrt(1+r.A2.*sin(th))','r','th')；

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

(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)

PagelOlExercise6

t=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(t,f)

PagelOlExercise7

xa=-l:0.l:l;ya=O:0.1:2;

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

z=x.*exp(x.A2+y.A2);

[zx,zy]=gradient(z,xa,ya);

f=sqrt(l+zx.A2+zy.A2);

s=0;

fori=2:length(xa)

forj=2:length(ya)

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

end

end

s

PagelOlExercise8

funl=inline('-(-x).A0.2.*cos(x),)；

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

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

PagelOlExercise9(以132为例)

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

h=O.1;x=0:h:32*pi;y=feval(fun,x);t1=trapz(x,y)

h=pi;x=O:h:32*pi;y=feval(fun,x);t2=trapz(x,y)%步长与周期一样，结果

失真

q1=quad(fun,0,32*pi)

q2=quadl(fun,0,32*pi)

PagelOlExercise10

⑵

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

functiond=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-dO)>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=inline(,x.A2*sin(x.A2-x-2)','x');

d=ex5_10_2f(fun,1.4,0.1,1e-3)

PagelOlExercise11

提示：f上升时，f'>O;f下降时，f'<示f极值,f'=0.

PagelOlExercise12

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

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_l2(inlineCl/sqrt(2*pi)*exp(-x.A2/2),),0,1,50)

PagelOlExercise13

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

quadl(fun,15,30)

PagelOlExercise14

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

PagelOlExercise15

利润函数fun=inline('(p-cO+k*log(M*exp(-a*p)))*M*exp(-a*p)',,p,);

求p使fun最大

PagelOlExercise16

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)%单位：千牛

PagelOlExercise17

clear;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%单位：百万元

PagelOlExercise18

曲面面积计算

Chapter6

Page121Exercise1

⑴

fun=inline('x+y','x','y')；

[t,y]=ode45(fun,[0123],1)%留意由于初值为y(O)=l,[O123]中。不

行缺

(3)

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

y⑴』y⑵,y⑵'=0.01*y⑵人2-2*y6+sin(t)

运行下列指令

clear;close;

fun=@(t,y)[y⑵;0.01*j^2)A2-2*y6+sin(t)];

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

plot(t,y(:,l))

(5)

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

y(l),=y(2),y(2)'=-mu*(y(l)A2-l)*y(2)-y(l)

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

clear;close;

fun=@(t,y,mu)[y⑵;-mu*(y⑴人2-l)*y⑵-y⑴];

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

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

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

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

Page121Excercise2

roots([l105413213750])

通解

Al*exp(-3*t)*cos(4*t)+A2*exp(-3*t)*sin(4*t)+A3*exp(-2*t)+A4*exp(-

t)+A5*t*exp(-t)

Pagel21Excercise3

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

*y(2)+0.5]7x'；y')；

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

所用节点很多

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

所用节点很少

Pagel21Excercise4

clear;

dfun=inline('[x(2);2*x(3)+x(1)-((l-l/82.45)*(x(l)+1/82.45))/(sqrt((x

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

+1-1/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)

AAAA',,

)3-(l/82.45*x(3))/(sqrt((x(l)+l-l/82.45)2+x(3)2))3]','t)x);

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

plot(x(:,l),x(:,3));

Page121Excercise5

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

clear;close;

fun=inline('2*x+yA2',"x,,'y,);

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

plot(x,y)

Pagel21Excercise6

clear;close;

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

[t,x]=ode45(fun,[010],0.

subplot(2,4,l)；plot(t,x)

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

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

[t,x]=ode45(fun,[010],0.1J],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.

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

[t,x]=ode45(fun,[0

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.

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

Page121Excercise7

微分方程T'=k(c-T),T(0)=20

dsolve(,DT=k*(c-T),,'T(0)=