數(shù)學(xué)實(shí)驗(yàn)課后習(xí)題解答

數(shù)學(xué)實(shí)驗(yàn)課后習(xí)題解答

配套教材：王向東戎海武文翰編著數(shù)學(xué)實(shí)驗(yàn)

王汝軍編寫

實(shí)驗(yàn)一曲線繪圖

【練習(xí)與思考】

畫出下列常見曲線的圖形。

以直角坐標(biāo)方程表示的曲線:

1.立方曲線＞=/

clear;

x=-2:0.1:2;

y=x.A3;

plot(x,y)

2.立方拋物線丁=1

clear;

y=-2:0.1:2;

x=y.A3;

plot(x,y)

gridon

3.高斯曲線ke

clear;

x=-3:0.1:3;

y=exp(-x.A2);

plot(x,y);

gridon

%axisequal

以參數(shù)方程表示的曲線

2

4.奈爾拋物線x=「，y=產(chǎn)曰=/)

clear;

t=-3:0.05:3;

x=t.A3;y=t.A2;

plot(x,y)

axisequal

gridon

5.半立方拋物線X=*,y=f3(y2=x3)

clear;

t=-3:0.05:3;

x=t.A2;y=t.A3;

plot(x,y)

%axisequal

gridon

3at3at2

6.迪卡爾曲線x=(尤3+y3―^axy-0)

1+產(chǎn)1+r

clear;

a=3;t=-6:0.1:6;

x=3*a*t./(l+t.A2);

y=3*a*t.A2./(l+t.A2);

plot(x,y)

233

7.蔓葉線%=%(y2=-^-)

1+，1+ECl—X

clear;

a=3;t=-6:0,1:6;

x=3*a*t.A2./(l+t.A2);

y=3*a*t.A3./(l+t.A2);

plot(x,y)

60廠-

40

20

0

-20

-40

-60

0123456789

8.擺線x=a(t-sint),y=b(l-cos/1)

clear;clc;

a=l;b=l;

t=0:pi/50:6*pi;

x=a*(t-sin(t));

y=b*(l-cos(t));

plot(x,y);

axisequal

gridon

31

9.內(nèi)擺線(星形線)x=6/costyy=asin1(x^+=a^)

clear;

a=l;

t=0:pi/50:2*pi;

x=a*cos(t).A3;

y=a*sin(t).A3;

plot(x,y)

1r-----------------------------------------------------------------------------------------------------------"J

0.8

0.6

0.4

0.2

0

-0.2

-0.4

-0.6

-0.8

-1-------------------------------------------------------1-------------------------------------------------------

?1-0.8-0.6-0.4-0.200.20.40.60.81

10.圓的漸伸線(漸開線)x=a(cosZ+rsint),y=a(sint-tcost)

clear;

a=l;

t=0:pi/50:6*pi;

x=a*(cos(t)+t.*sin(t));

y=a*(sin(t)+t.*cos(t));

plot(x,y)

gridon

11.空間螺線x=acost.y=bsint,z=ct

clear

a=3;b=2;c=l;

t=0:pi/50:6*pi;

x=a*cos(t);

y=b*sin(t);

z=c*t;

plot3(x,y,z)

gridon

以極坐標(biāo)方程表示的曲線：

12.阿基米德線r—a(p,r>0

clear;

a=l;

phy=0:pi/50:6*pi;

rho=a*phy;

polar(phy,rho/r-*r)

13.對數(shù)螺線r=e"。

clear;

a=0.1;

phy=0:pi/50:6*pi;

rho=exp(a*phy);

polar(phy,rho)

14.雙紐線r2=a2cos2(p((x2+y2)2=a2(x2-y2))

clear;

a=l;

phy=-pi/4:pi/50:pi/4;

rho=a*sqrt(cos(2*phy));

polar(phy,rho)

holdon

polar(phy<rho)

15.雙紐線r~=a~sin2(p((x2+y2)2=2a2xy)

clear;

a=l;

phy=0:pi/50:pi/2;

rho=a*sqrt(sin(2*phy));

polar(phy,rho)

holdon

polar(phy<rho)

16.四葉玫瑰線r=6/sin2(p,r>0

clear;close

a=l;

phy=0:pi/50:2*pi;

rho=a*sin(2*phy);

polar(phy,rho)

17?二葉玫瑰線r—asin3(p,r>0

clear;close

a=l;

phy=0:pi/50:2*pi;

rho=a*sin(3*phy);

polar(phy,rho)

18.三葉玫瑰線r-tzcos3^,r>0

clear;close

a=l;

phy=0:pi/50:2*pi;

rho=a*cos(3*phy);

polar(phy,rho)

實(shí)驗(yàn)二極限與導(dǎo)數(shù)

【練習(xí)與思考】

1.求下列各極限

(1)lim(l--)n(2)limVn3+3"(3)+2-27^71+4)

M—＞00幾〃一＞8〃一＞8

clear;

symsn

yl=limit((l-l/n)An,n,inf)

y2=limit((nA3+3An)A(l/n),n,inf)

y3=limit(sqrt(n+2)-2*sqrt(n+l)+sqrt(n),n,inf)

yl=l/exp(1)

y2=3

y3=0

(4)lim(-------------)(5)limxcot2x(6)lim(Jx2+3x-x)

XiX~~1X-lXTOXT8

clear;

symsx;

y4=limit(2/(xA2-l)-l/(x-l),x；l)

y5=limit(x*cot(2*x),x,0)

y6=limit(sqrt(xA2+3*x)-x,x,inf)

y4=-1/2

y5=1/2

y6=3/2

(7)lim(cos—)J(8)lim(-———)(9)lim^1+—1

XT8XXTlXeX-1I。X

clear;

symsxm

y7=limit(cos(m/x),x,inf)

y8=limit(l/x-l/(exp(x)-l),x,l)

y9=limit(((l+x)A(l/3).l)/x,x,0)

y7=1

y8=(exp(1)-2)/(exp(1)-1)

y9=1/3

2.考慮函數(shù)

/(x)=3x2sin(x3),-2<x<2

作出圖形，并說出大致單調(diào)區(qū)間；使用diff求廣(x),并求/(x)確切的單調(diào)區(qū)間。

clear;close;

symsx;

f=3*xA2*sin(xA3);

ezplot(邙?2,2])

gridon

大致的單調(diào)增區(qū)間：[2?1刃,[?1.3,1.2],[1.7,2];

大致的單點(diǎn)減區(qū)間：[-1.7,-L3],[L2,L7];

fl=diff(f,xj)

ezplot(fl9[-2,2])

line([-5,5],[0,0])

gridon

axis([-2.1,2.l,?60,120])

fl=

6*x*sin(xA3)+9*xA4*cos(xA3)

用fzero函數(shù)找f\x)的零點(diǎn),即原函數(shù)/(x)的駐點(diǎn)

xl=fzero(,6*xH：sin(xA3)+9*xA4*cos(xA3)\[-2,-1.7])

x2=fzero(/6*x*sin(xA3)+9*xA4*cos(xA3)\[-1.7,-1.5])

x3=fzero(*6*x*sin(xA3)+9*xA4*cos(xA3)\[-1.5<l.l])

x4=fzero(/6*x*sin(xA3)+9*xA4*cos(xA3)\0)

x5=fzero(*6*x*sin(xA3)+9*xA4*cos(xA3)\[l,l-5])

x6=fzero(*6*x*sin(xA3)+9*xA4*cos(xA3)\[1.5,1.7])

x7=fzero(*6*x*sin(xA3)+9*xA4*cos(xA3)\[1.7,2])

xl=

-1.9948

x2=

-1.6926

x3=

-1.2401

x4=

0

x5=

1.2401

x6=

1.6926

x7=

1.9948

確切的單調(diào)增區(qū)間：[-1.9948,-1.6926],[-1.2401,1.2401],[1.6926,1.9948]

確切的單調(diào)減區(qū)間：[-2,-1.9948],[-1.6926,-1.24011,[1.2401,1.6926],[1.9948,2]

3.對于下列函數(shù)完成下列工作，并寫出總結(jié)報(bào)告，評論極值與導(dǎo)數(shù)的關(guān)系，

⑴作出圖形，觀測所有的局部極大、局部極小和全局最大、全局最小值點(diǎn)的粗略

位置；

(il)求_f(x)所有零點(diǎn)(即/(x)的駐點(diǎn))；

(iii)求出駐點(diǎn)處/(x)的二階導(dǎo)數(shù)值；

(iv)用fmin求各極值點(diǎn)的確切位置；

(v)局部極值點(diǎn)與/'(x)J"(x)有何關(guān)系？

(1)/(%)=x2sin(x2-x-2),xe[-2,2]

(2)/(x)=3x5-20x3+1(),XG[-3,3]

(3)/(x)=|x3-x2-x-2|,xe[0,3]

clear;close;

symsx;

f=xA2*sin(xA2-x-2)

ezplot(f,[-2,2])

gridon

f=

x人2*sin(x人2-x-2)

局部極大值點(diǎn)為：-1.6,局部極小值點(diǎn)為為：-0.75,-1.6

全局最大值點(diǎn)為為：-1.6,全局最小值點(diǎn)為：-3

fl=diff(f,x,l)

ezplot(fl,[-2,2])

line([-5,5],[0,0])

gridon

axis([-2,1,2.1,-6,201)

fl=

2*x*sin(xA2-x-2)+xA2*cos(x^2-x-2)*(2*x-1)

用fzero函數(shù)找f(x)的零點(diǎn)，即原函數(shù)/(x)的駐點(diǎn)

xl=fzero(/2*x*sin(xA2-x-2)+xA2*cos(xA2-x-2)*(2*x-l)\[-2,-1.2])

x2=fzero(,2*x*sin(xA2-x-2)+xA2*cos(xA2-x-2)*(2*x-l)\[-1.2,-0.5])

x3=fzero(/2*x*sin(xA2-x-2)+xA2*cos(xA2-x-2)*(2*x-l),,[-0.5,1.2])

x4=fzero(,2*x*sin(xA2-x-2)+xA2*cos(xA2-x-2)*(2*x-l)\[1.2,2])

xl=

-1.5326

x2=

-0.7315

x3=

-3.2754e-027

x4=

1.5951

ff=@(x)X.A2.*sin(x.A2-x-2)

ff(-2),ff(xl),ff(x2),ff(x3),ff(x4),ff⑵

ff=

@(x)x.A2.*sin(x.A2-x-2)

ans=

-3.0272

ans=

2.2364

ans=

-0.3582

ans

-9.7549e-054

ans=

-2.2080

ans=

0

實(shí)驗(yàn)三級數(shù)

【練習(xí)與思考】

1.用taylor命令觀測函數(shù)y=/(x)的Maclaurin展開式的前幾項(xiàng),然后在同一坐

標(biāo)系里作出函數(shù)y=/(x)和它的Taylor展開式的前幾項(xiàng)構(gòu)成的多項(xiàng)式函數(shù)的圖形，

觀測這些多項(xiàng)式函數(shù)的圖形向y=/(x)的圖形的逼近的情況

(1)f(x)=arcsinx

clear;

symsx

y=asin(x);

yl=taylor(y,0,l)

y2=taylor(y,0,5)

y3=taylor(y,0,10)

y4=taylor(y,0,15)

y=subs(y,x);

yl=subs(yl,x);

y2=subs(y2,x);

y3=subs(y3,x);

y4=subs(y4,x);

plot(x,y,x,yl,':,,x,y2,'-.',x,y3,,-',x,y4,':','linewidth',3)

yl=

0

y2=

x-3/6+x

y3=

AAA

(35*x9)/1152+(5*x7)/112+(3*x、5)/40+x3/6+x

y4=

(231*xA13)/13312+(63*xAll)/2816+(35*x"9)/1152+(5*xA7)/112+

A

(3*x5)/40+x人3/6+x

(2)f(x)=arctanx

clear;

symsx

y=atan(x);yl=taylor(y,0,3)

y2=taylor(y,0,5),y3=taylor(y,0,10),y4=taylor(y,0,15)

x=-l:0.1:l;

y=subs(y,x);yl=subs(yl,x);y2=subs(y2,x);

y3=subs(y3,x);y4=subs(y4,x);

plot(x,y,x,yl,':',x,y2,'-.',x,y3,'-,,x,y4,':','linewidth',3)

yi=

X

y2=

x-xA3/3

y3=

x-9/9-xA7/7+xA5/5-xA3/3+x

y4=

AAA

x人13/13-xll/ll+x人9/9-x7/7+x5/5-x人3/3+x

4

⑶/(x)=e，

clear;

symsx

y=exp(xA2);

yl=taylor(y,0,3)

y2=taylor(y,0,5)

y3=taylor(y,0,10)

y4=taylor(y,0,15)

x=-l:0.1:l;

y=subs(y,x);

yl=subs(yl,x);

y2=subs(y2,x);

y3=subs(y3,x);

y4=subs(y4,x);

plot(x,y,x,yl,3,'linewidth',3)

yi=

x-2+1

y2=

x八4/2+x-2+1

y3=

xA8/24+xA6/6+xA4/2+xA2+1

y4=

AAAAA

x人14/5040+x12/720+x10/120+x8/24+x6/6+x人4/2+x2+1

(4)/(x)=sin2x

clear;

symsx

y=sin(x)A2;

yl=taylor(y,0,l)

y2=taylor(y,0,5)

y3=taylor(y,0,10)

y4=taylor(y,0,15)

x=-pi:0.1:pi;

y=subs(y,x);

yl=subs(yl,x);

y2=subs(y2,x);

y3=subs(y3,x);

y4=subs(y4,x);

plot(x,y,x,yl,3,'linewidth',3)

yi=

o

y2=

x-2-X-4/3

y3

-xA8/315+(2*xA6)/45-xA4/3+xA2

y4=

(4*xA14)/42567525-(2*xA12)/467775+(2*xA10)/14175-xA8/315+

AA

(2*x6)/45-x人4/3+x2

5

(5)/(x)=--

X

clear;

symsx

y=exp(x)/(l-x);

yl=taylor(y,0,3)

y2=taylor(y,0,5)

y3=taylor(y,0,10)

y4=taylor(y,0,15)

x=-l:0.1:0;

y=subs(y,x);

yl=subs(yl,x);

y2=subs(y2,x);

y3=subs(y3,x);

y4=subs(y4,x);

plot(x,y,x,yl,':',x,y2,'-.',x,y3,'linewidth',3)

yi=

(5*xA2)/2+2*x+1

y2=

(65*xA4)/24+(8*xA3)/3+(5*xA2)/2+2*x+1

y3=

(98641*x7)/36288+(109601*xA8)/40320+(685*xA7)/252+

(1957*x-6)/720+(163*xA5)/60+(65*xA4)/24+(8*xA3)/3+(5*xA2)/2

+2*x+1

y4=

(47395032961*xA14)(8463398743*xA13)/3113510400+

(260412269*xA12)/95800320+(13563139*xAll)/4989600+

AA

(9864101*x10)/3628800+(98641*x人9)/36288+(109601*x8)/40320+

(685*x”)/252+(1957*xA6)/720+(163*xA5)/60+(65*x7)/24+

(8*xA3)/3+(5*xA2)/2+2*x+1

(6)/(x)=ln(x+71+x2)

clear;

symsx

y=log(x+sqrt(l+xA2));

yl=taylor(y,0,3)

y2=taylor(y,0,5)

y3=taylor(y,0,10)

y4=taylor(y,0,15)

x=-l:0.1:l;

y=subs(y,x);

yl=subs(yl,x);

y2=subs(y2,x);

y3=subs(y3,x);

y4=subs(y4,x);

plot(x,y,x,yl,':',x,y2,'-.',x,y3,'-',x,y4,':','linewidth',3)

yi=

X

y2=

x-xA3/6

y3=

(35*xA9)/1152-(5*xA7)/112+(3*xA5)/40-xA3/6+x

y4=

(231*xA13)/13312-(63*xAll)/2816+(35*xA9)/1152-(5*xA7)/112+

A

(3*x5)/40-x八3/6+x

4

oo1"2k

2.求公式=—#=12?)中的數(shù)恤#=4,5,6,7,8的值.

n=i〃叫

k=[45678];

symsn

symsum(l./n?八(2*k),l,inf)

ans=

A

[pi人8/9450,pi人10/93555,(691*pi12)/638512875,(2*pi人14)/18243225,

(3617*piA16)/325641566250]

81

3.利用公式=e來計(jì)算e的近似值。精確到小數(shù)點(diǎn)后100位,這時(shí)應(yīng)計(jì)算到

n=04

這個(gè)無窮級數(shù)的前多少項(xiàng)?請說明你的理由.

解：Matlab代碼為

clear;clc;close

epsl=1.0e-100;

ep=l;fn=l;a=l;n=l;

whileep>epsl

a=a+fn;

n=n+l;

fn=fn/n;

ep=fn;

end

fn

vpa(a,100)

n

8.3482e-101

ans=

2.71828182845904553488480814849026501178741455078125

n=

70

精確到小數(shù)點(diǎn)后100位，這時(shí)應(yīng)計(jì)算到這個(gè)無窮級數(shù)的前71項(xiàng)，理由是誤差小于

1()的負(fù)10()次方，需要最后一項(xiàng)小于10的負(fù)1()0次方，由上述循環(huán)知n=7()時(shí)最

后一項(xiàng)小于10的負(fù)100次方，故應(yīng)計(jì)算到這個(gè)無窮級數(shù)的前71項(xiàng).

4.用練習(xí)3中所用觀測法判斷下列級數(shù)的斂散性

clear;clc;

epsl=0.000001;

N=50000;p=1000;

symsn

Un=l/(nA2+nA3);

sl=symsum(Un,l,N);

s2=symsum(Un,l,N+p);

sa=vpa(s2-sl);

sa=setstr(sa);

sa=str2num(sa);

fprintfC級數(shù)，)

disp(Un)

ifsa<epsl

dispC收斂，)

else

dispC發(fā)散)

end

級數(shù)1/(nA3+n人2)收斂

clear;close

symsn

s=[]；

fork=l:l()0

s(k)=symsum(l/(nA3+nA2),l,k);

end

plot(s/.*)

0.66

0.64

0.62

0.6

0.58

0.56

0.54

0.52

0.5

0102030405060708090100

clear;clc;

epsl=O.OOOOOl;

N=50000;p=1000;

symsn

Un=l/(n*2An);

sl=symsum(UnJ,N);

s2=symsum(Un,l,N+p);

sa=vpa(s2-sl);

sa=setstr(sa);

sa=str2num(sa);

fprintfC級數(shù)？

disp(Un)

ifsa<epsl

dispC收斂，)

else

dispC發(fā)散，)

end

級數(shù)1/(2人n*n)收斂

clear;close

symsn

s=[]；

fork=l:l()()

s(k)=symsum(l/(2An*n),l,k);

end

plot(s；J)

0.7

0.68

0.66

0.64

0.62

0.6

0.58

0.56

0.54

0.52

0102030405060708090100

8I

(3)ysin-

〃=]n

clear;clc;

epsl=O.OOOOOOOOOOOOOl;

N=50000;p=100;

symsn

Un=l/sin(n);

sl=symsum(Un,l,N);

s2=symsum(Un,l?N+p);

sa=vpa(s2-sl);

sa=setstr(sa);

sa=str2num(sa);

fprintfC級數(shù)，

disp(Un)

ifabs(sa)<epsl

dispC收斂，)

else

disp('W)

end

級數(shù)l/sin(n)發(fā)散

clear;close

symsn

s=[]；

fork=l:l()0

s(k)=symsum(l/sin(n)J,k);

end

plot(s，')

發(fā)散

⑷落

n=ln

clear;clc;

epsl=0.0()00001;

N=50000;p=1000;

symsn

Un=log(n)/(nA3);

sl=symsum(Un,l,N);

s2=symsum(Un,l?N+p);

sa=vpa(s2-sl);

sa=setstr(sa);

sa=str2num(sa);

fprintfC級數(shù)，

disp(Un)

ifsa<epsl

dispC收斂，)

else

disp('W)

end

級數(shù)log(n)/nA3收斂

clear;close

symsn

s=[]；

fork=l:l()0

s(k)=symsum(log(n)/nA3,l,k);

end

plot(s，')

⑸*

n=\n

clear;close

symsn

s=[];he=0;

fork=l:100

he=he+factorial(k)/kAk;

s(k)=he;

end

plot(s/J)

1.9

1.8

1.7

1.6

1.5?

1.4

1.3

1.2

1.1

0102030405060708090100

1

(6)y--------

clear;clc;

epsI=0.0000001;

N=50000;p=1000;

symsn

Un=l/log(n)An;

sl=symsum(Un,35N);

s2=symsum(Un,3^N+p);

sa=vpa(s2-sl);

sa=setstr(sa);

sa=str2num(sa);

fprintfC級數(shù)，)

disp(Un)

ifsa<epsl

dispC收斂)

else

dispC發(fā)散)

end

級數(shù)l/log(n)An收斂

clear;close

symsn

s=[]；

fork=3:100

s(k)=symsum(l/log(n)An,3?k);

end

plot(s,'.')

1.4[

1.2

1

0.8

?

0.6

0.4

0.2

0.

0102030405060708090100

⑺￡1

n=\Inn

clear;clc;

epsl=O.OOOOOOl;

N=50000;p=100;

symsn

Un=l/(log(n)*n);

sl=symsum(Un,3?N);

s2=symsum(Un,3?N+p);

sa=vpa(s2-sl);

sa=setstr(sa);

sa=str2num(sa);

fprintf。級數(shù)，

disp(Un)

if(sa)<epsl

dispC收斂)

else

dispC發(fā)散)

end

級數(shù)1/(n*log(n))發(fā)散

clear;close

symsn

s=[]；

fork=3:300

s(k)=symsum(l/(n*log(n)),2,k);

end

⑻6署

w=ln十I

clear;clc;

epsl=0.00000()l;

N=50000;p=100;

symsn

Un=(-l)An*n/(nA2+l);

sl=symsum(Un,3?N);

s2=symsum(Un,3,N+p);

sa=vpa(s2-sl);

sa=setstr(sa);

sa=str2num(sa);

fprintfC級數(shù)，

disp(Un)

if(sa)<epsl

dispC收斂，)

else

dispC發(fā)散，)

end

級數(shù)((-1)人n*n)/(rT2+1)收斂

clear;close

symsn

s=[]；

fork=3:300

s(k)=symsum((-l)An*n/(nA2+l),2,k);

end

plot(s；/)

實(shí)驗(yàn)四積分

【練習(xí)與思考】

1.（不定積分）用int計(jì)算下列不定積分,并用diff驗(yàn)證

2rdx3

fxsinxtZx,[—————Je'+lIarcsinxdx,[secxdx

JJ1+cosx

解：Matlab代碼為：

symsx

yl=x*sin(xA2);

y2=l/(l+cos(x));

y3=l/(exp(x)+l);

y4=asin(x);

y5=sec(x)A3;

fl=int(yl)

f2=int(y2)

f3=int(y3)

f4=int(y4)

f5=int(y5)

dy=simplify(diff([fl;f2;f3;f4;f5]))

dy=

x*sin(xA2)

tan(x/2)A2/2+1/2

1/(exp(x)+1)

asin(x)

(cot(pi/4+x/2)*(tan(pi/4+x/2)人2/2+1/2))/2+1/(2*cos(x))+

tan(x)A2/cos(x)

fl=

-cos(x人2)/2

f2=

tan(x/2)

f3=

x-log(exp(x)+1)

f4=

x*asin(x)+(1-x人2)人(1/2)

f5=

log(tan(pi/4+x/2))/2+tan(x)/(2*cos(x))

2.(定積分)用trapz,quad,int計(jì)算下列定積分

，fx'tZx,exsin(2x)Jx,e~x2dx

Jo%JoJoJo

解：Matlab代碼為

clear;

x=(0+eps):0.05:l;

yl=sin(x)Jx;

fl=trapz(x,yl)

fl=0.9460

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

fl2=quad(funl,0+eps,l)

fl2=0.9461

fl3Kpa(int('sin(x)/x',0,l),5)

fl3=0.94608

r2v2

3.(橢圓的周長)用定積分的方法計(jì)算橢圓三+'=1的周長

解：橢圓的參數(shù)方程為['="os’

[y=2sinr

由參數(shù)曲線的弧長公式得

s=/+y\t)2dt=『>/9sin2/+4cos2tdt=『>/5sin2t+4dt

Matlab代碼為

s=vpa(int(,sqrt(5*sin(t)A2+4)*/t\0,2*pi),5)

s=

15.865

4.(二重積分)計(jì)算數(shù)值積分\\(\+x+y)dxdy

x2+y2<,2y

解：fxy=@(x,y)l+x+y;ylow=@(x)l-sqrt(l-x.A2);yup=@(x)l+sqrt(l-x.A2);

s=quad2d(fxy<l,l,ylow,yup)

s=

6.2832

或符號積分法：

symsxy

xi=int(l+x+y,y,l-sqrt(l-xA2),l+sqrt(l-xA2));

s=int(xi,x,-l,l)

s=

2*pi

5.(假奇異積分)用trapz,quad8計(jì)算積分1,"cosxcZr,會(huì)出現(xiàn)什么問題?

分析原因，并求出正確的解。

解：Matlab代碼為

clear

x=-l:0.05:l;

y=x.A(l/3).*cos(x);

sl=trapz(x,y)

fun5=@(x)x.A(l/3).*cos(x);

s2=quad(fun5,-l,l)

int(,xA(l/3)*cos(x),/x\-l,l)

si=

0.9036+0.5217i

s2=

0.9114+0.5262i

Warning:Explicitintegralcouldnotbefound.

ans=

int(xA(1/3)*cos(x),x=-1..1)，原函數(shù)不存在，不能用int函數(shù)運(yùn)算。

用梯形法和辛普森法計(jì)算數(shù)值積分時(shí)，由于對負(fù)數(shù)的開三次方運(yùn)算結(jié)果為復(fù)

數(shù)，所以導(dǎo)致結(jié)果錯(cuò)誤且為復(fù)數(shù)；

顯然被積函數(shù)為奇函數(shù)，在對稱區(qū)間上的積分等于0,此時(shí)可以這樣處理：

(1)重新定義被積函數(shù)

%fun5.m

functiony=fun5(x)

[m,n]=size(x);

fork=l:m

forl=l:n

y(k,l)=nthroot(x(k,l),3)*cos(x(kj));

end

end

end

用辛普森法：

s=quad(,fun5\-l,l)

s=

0

用梯形法

clear;

x=-l:0.01:l;

y=fun5(x);

s=trapz(x,y)

s=

-1.3878e-017

6.(假收斂現(xiàn)象)考慮積分/(Z)=(：sinx團(tuán)，

(1)用解析法求/伙)；

clear;

symsxk;

Ik=int(abs(sin(x))909k*pi)

Warning:Explicitintegralcouldnotbefound.

Ik=

int(abs(sin(x))rx=0.,pi*k)

(2)分別用trapz,quad和quad8求/(4),/(6)和/(8),發(fā)現(xiàn)什么問題？

clear;

fork=4:2:8;

x=0:pi/1000:k*pi;

y=abs(sin(x));

trapz(x,y)

end

ans

8.0000

ans=

12.0000

ans=

16.0000

fork=4:2:8

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

quad(fun6,0,k*pi)

end

ans=

8.0000

ans=

12.0000

ans=

16.0000

7.(Simpson積分法)編制一個(gè)定步長Simpson法數(shù)值積分程序.計(jì)算公式為

h

/=s”=](力+4力+2力+4力+…+2/?,,+4/?+fn+i)

其中〃為偶數(shù)，力=-―-fi=f(a+(z-l)A),z=1,2,+1.

ny

解：Matlab代碼為

%fun7.m

functiony=fun7(Lname,a,b,n)

%f_name為被積函數(shù)

%[a,b]為積分區(qū)間

%n為偶數(shù)，用來確定步長h=(b-a)/n

ifmod(n,2)~=0

disp(rn必須為偶數(shù)')

return;

end

ifnargin<4

n=100;

end

ifnargin<3

dispC請輸入積分區(qū)間)

end

ifnargin==0

dispC^rror1)

end

h=(b-a)/n;

x=a:h:b;

s=0;

fork=l:n+l

ifk==l||k==(n+l)

xishu=l;

elseifmod(k,2)==0

xishu=4;

else

xishu=2;

end

s=s+feval(fLname,x(k))*xishu;

end

y=s*h/3;

end

8.(廣義積分)計(jì)算廣義積分

exp(-x2).ritan(x)?Isinx,

；

-----4—dx,―^—dx,,ax

L1+xJ。6\-x2

并驗(yàn)證公式

x2

exp(-----)

sinx

dx=-

J-8J2萬x2

解：Matlab代碼為

clear;

symsx

sl=vpa(int(exp(-xA2)/(l+xA4)<inf,inf),5)

s2=quad(@(x)tan(x)./sqrt(x),0,l)

s3=quad(@(x)sin(x)./sqrt(l-x.A2),0,l)

A

s4=vpa(int(exp(-x2/2)/sqrt(2*pi),-inf,inf),5)

s5=int(sin(x)./x,0+eps,inf)

si=

1.4348

s2=

0.7968

s3=

0.8933

s4=

1.0

s5=

pi/2-sinint(1/4503599627370496)

實(shí)驗(yàn)五二元函數(shù)的圖形

【練習(xí)與思考】

1.畫出空間曲線+在_30<%”30范圍內(nèi)的圖形,并畫

出相應(yīng)的等高線。

clear;

x=-30:0.5:30;y=-30:0.5:30;

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

Z=10*sin(sqrt(X.A2+Y.A2))./sqrt(l+X.A2+Y.A2);

mesh(X,Y,Z)

2.根據(jù)給定的參數(shù)方程，繪制下列曲面的圖形。

a)橢球面x=3coszzsinv,y=2coswcosv,z=sinw

clear;

u=0:pi/50:2*pi;

v=0:pi/50:pi;

[U,V]=meshgrid(u,v);

x=3*cos(U).*sin(V);

y=2*cos(U).*cos(V);

z=sin(U);

mesh(x,y,z)

b)橢圓拋物面x=3〃sin匕y=2wcosv,z=4u2

clear;

u=0:pi/50:pi/4;

v=0:pi/50:2*pi;

[U,V]=meshgrid(u,v);

x=3*U.*sin(V);

y=2*U.*cos(V);

z=4*U.A2;

mesh(x,y,z)

axisequal

C)單葉雙曲面無二3sec〃sin匕y=2sec〃cosu,z4tanu

clear;

u=0:pi/15:pi;

v=0:pi/15:2*pi;

[U,V]=meshgrid(u,v);

x=3*sec(U).*sin(V);

y=2*sec(U).*cos(V);

z=4*tan(U);

mesh(x,y,z)

-20-40

22

d)雙曲拋物面x=u,y=v,z="§"

clear

u=-3:0.1:3;

[U,V]=meshgrid(u);

x=U;

y=v；

z=(U.A2-V.A2)/3;

mesh(x9y9z)

e)旋轉(zhuǎn)面x=In〃sin匕y=In〃cosv,z=u

clear;

u=l:0.1:5;

v=0:pi/30:2*pi;

[U,V]=meshgrid(u,v);

x=log(U).*sin(V);

y=log(U).*cos(V);

z=U;

mesh(x,y,z)

axisequal

f)圓錐面x=〃sin匕y=ucosv,z=u

clear;

u=-5:0.1:5;

v=0:pi/30:2*pi;

[U,V]=meshgrid(u,v);

x=(U).*sin(V);

y=(U).*cos(V);

z=U;

mesh(x,y,z)

axisequal

g)環(huán)面x=(3+0.4cosu)cos匕y=(3+0.4cos〃)sin匕z=0.4sinv

clear;

u=0:pi/30:2*pi;

v=u;

[U,V]=meshgrid(u,v);

x=(3+0.4*cos(U)).*cos(V);

y=(3+0.4*cos(U)).*sin(V);

z=0.4*sin(V);

mesh(x9y9z)

h)正螺面x=wsinv,y=ucosv,z=4v

clear;

u=0:pi/30:pi;

v=0:pi/30:10*pi;

[U,V]=meshgrid(u,v);

x=U.*sin(V);

y=U.*cos(V);

z=4*V;

mesh(x,y,z)

colorbar

3.在一丘陵地帶測量搞程,x和y方向每隔100米測一個(gè)點(diǎn),得高程

見表5-2,試擬合一曲面，確定合適的模型，并由此找出最高點(diǎn)和該點(diǎn)的高

程.

表5-2高程數(shù)據(jù)

100200300400

100636697624478

200698712630478

300680674598412

400662626552334

clc;clear;

xl=[100100100100200200200200300300300300400400400400];

x2=[100200300400100200300400100200300400100200300400];

y=[636698680662697712674626624630598552478478412334]';

x=[xl\x2r];

x0=[l1111];

beta=lsqcurvefit(*heigh\xO,x,y)

%繪圖：

al=100:5:400;

a2=al;

[xxl,xx2]=meshgrid(al,a2);

Z=beta(l)+beta(2)*xxl+beta(3)*xx2+beta(4)*xxl.A2+beta(5)*xx2.A2;

mesh(xxl,xx2,Z)

Localminimumpossible.

Isqcurvefitstoppedbecausethefinalchangeinthesumofsquares

relativeto

itsinitialvalueislessthanthedefaultvalueofthefunction

tolerance.

beta=

Columns1through5

538.43751.49010.6189-0.0046-0.0017

contour(xxl,xx2,Z,30),colorbar

%計(jì)算最高點(diǎn)及高程

x0=[100,100];

options=optimset('largescale','off');

%設(shè)置下界

lb=[0,0];

%無上界

ub=[];

[x,fval]=fmincon(,height\xO,[],[],[],[],lb,ub,[],options)

Warning:OptionsLargeScale=1off'andAlgorithm=*trust-region-

reflective*conflict.

IgnoringAlgorithmandrunningactive-setalgorithm.Toruntrust-

region-reflective,set

LargeScale='on'.Torunactive-setwithoutthiswarning,use

Algorithm=factive-set,.

>Infminconat445

Localminimumpossible.Constraintssatisfied.

fminconstoppedbecausethepredictedchangeintheobjective

function

islessthanthedefaultvalueofthefunctiontoleranceand

constraints

weresatisfiedtowithinthedefaultvalueoftheconstraint

tolerance.

Noactiveinequalities.

x=

161.9676182.0320

fval=

-715.4403

400.700

350I650

300600

550

250

500

200

U450

150

Ho

100

100150200250300350400

heigh和height兩個(gè)函數(shù)分別定義如下：(應(yīng)寫在m文件中)

%heigh.m

functionf=heigh(beta,xdata)

xxl=xdata(:,1);

xx2=xdata(:,2);

f=beta(1)+beta(2)*xxl+beta(3)*xx2+beta(4)*xxl.八2+beta(5)*xx2.*xxl+beta(6

)*xx2.人2;

end

%height.m

functiony=height(x)

y=-(434.0000+1.9079*x(l)+1.0366*x(2)-0.0017*x(1).A2-0.0046*x(2).*x(l)-

0.0017*x(2).A2);

end

實(shí)驗(yàn)六多元函數(shù)的極值

【練習(xí)與思考】

1.求z=/+y4-4q+1的極值，并對圖形進(jìn)行觀測。

解：Maltab代碼為

symsxy;

z=xA4+yA4-4*x*y+l;

dzx=diff(z,x);

dzy=diff(z,y);

[x,y]=solve(dzx,dzy,x,y)

X=

0

1

-1

(-1)人(3/4)

一(一1)人(3/4)

-(-1)A(3/4)*i

(-1)人(3/4)*i

y=

0

1

-1

(-l)A(l/4)

-(-l)A(l/4)

i

-i

(-1)A(1/4)*i

-(-1)A(1/4)*i

經(jīng)計(jì)算可知，函數(shù)的駐點(diǎn)為(0,0)、(1,1)、(-1,-1)

ezmeshc(z,[-2,2<2,2])

從圖形上觀測可知，(1』)、(?1,?1)為極值點(diǎn)，(0,0)不是極值點(diǎn)。

clear

symsxy;

z=xA4+yA4-4*x*y+l;

dzx=diff(z,x);

A=diff(z,x,2)

B=diff(dzx,y)

C=diff(z,y,2)

A=

12*xA2

B=

-4

C=

12*yA2

由判別法可知(1,1)、(-1,-1)均為極小值點(diǎn)。

2.求函數(shù)在圓周/+y2=]的最大值和最小值。

解：構(gòu)造Lagrange函數(shù)

L(x,y)=/+2y2+2(x2+y2-1)

求Lagrange函數(shù)的自由極值.先求L關(guān)于x,y,/l的一階偏導(dǎo)數(shù)，再解正規(guī)方程可得

所求的極值點(diǎn)，Matlab代碼為

clear;

symsxyk

L=xA2+2*yA2+k*(xA2+yA2-l);

dlx=diff(L,x);

dly=diff(L,y);

dlk=diff(L,k);

s=solve(dlx,dly,dlk,x,y,k);

k=sk

x=s.x*

y=s?y'

k=

[T,-2,-1,-2]

x=

[1,0,-1/0]

y=

[o,i,o,-i]

t=0:pi/50:2*pi;

x=cos(t);

y=sin(t);

z=x.A2+2*y.A2;

plot3(x,y,z)

2

1.8

1.6

1.4

1.2

1

1

0.51

00-5

05-0.5

-1-1

可得點(diǎn)(1,0)、(0,1)(-1,0).(0,-1)為函數(shù)的條件極值點(diǎn)，經(jīng)判斷函

數(shù)/(%>)=/+2/在(1,0)、(-1,0)取得極小值，在(0,1)、(0,-1)取

得極大值。

3.在球面V+y2+z2=l求出與點(diǎn)(3,1,4)距離最近和最遠(yuǎn)點(diǎn)。

解：設(shè)球面上的點(diǎn)為(x,y,z),則此點(diǎn)與點(diǎn)(3,1,-1)的距離為

d(x,y,z)=(x-3)2+(y-1)2+(z+l)2且(x,y,z)滿x2+y2+z2=\；構(gòu)造

Lagrange函數(shù)

L(x,y,z,A)=(x—3)2