中南大学数学建模实验作业

数学实验与数学建模实验报告

学 院： 航空航天学院

专业班级： 探控1201班 姓 名： 安前松 学 号： 4201120124

完成时间： 2014 年1 月6日

1

承 诺 书

本人承诺所呈交的数学实验与数学建模作业都是本人通过学习自行进行编程独立完成，所有结果都通过上机验证，无转载或抄袭他人，也未经他人转载或抄袭。若承诺不实，本人愿意承担一切责任。

承诺人：

2014年1月 6日

注意事项如下：

1、上机时间：第11周到第18周星期六晚上：7：00-8：30； 2、上机地点：新校区数学与统计学院二楼数学实验室；

3、交作业时间：第19周的星期五（2014年1月 6日）交到新校区数学与统计学院二楼数学实验室数学实验室办公室；

4、报告所有结果交电子打印稿到新校区数学与统计学院二楼数学实验室办公室，并将电子文档发送到邮箱：xuanyunqin@163.com（word文档命名：姓名＋学号＋数学实验作业）

数学实验学习体会

（每个人必须要写1500字以上，占总成绩的20%）

2

实验一 图形的画法

1. 做出下列函数的图像：

22(1)y(x)?xsin(x?x?2)，?2?x?2（分别用plot、fplot） >>

x=-2:0.01:2;

y=x.^2.*sin(x.^2-x-2);

plot(x,y)

（2）x/9?y/25?1（用参数方程） >>

t=0:pi/50:2*pi; x=3*cos(t);y=5*sin(t);

plot(x,y)

22

(3) 在同一图形窗口中，画出四幅不同图形（用subplot命令）：

y1?cos(x)，y2?sin(x?pi/2)，y3?x2cos(x?pi)，y4?esin(x)（x?[0,2?]）

>>

x=0:pi/50:2*pi; subplot(2,2,1) y1=cos(x); plot(x,y1); subplot(2,2,2) y2=sin(x-pi/2); plot(x,y2); subplot(2,2,3)

y3=x.^2.*cos(x-pi); plot(x,y3); subplot(2,2,4) y4=exp(sin(x)); plot(x,y4)

3

2 作出极坐标方程为r?2(1?cost)的曲线的图形. >>

t=0:pi/50:2*pi; r=2*(1-cos(t)); polar(r,t)

3 作出极坐标方程为r?et/10的对数螺线的图形. >>

t=0:pi/50:2*pi; r=exp(t/10); x=r.*cos(t); y=r.*sin(t); plot(x,y)

?x?4cost,?4 绘制螺旋线?y?4sint,在区间［0，4?］上的图形.在上实验中，显示坐标

?z?t?轴名称。 >>

t=0:pi/50:4*pi; x=4*cos(t); y=4*sin(t); z=t;

plot3(x,y,z)

5 作出函数z??xye?x?y的图形.

22

4

>>

x=-1:0.01:1;

[x,y]=meshgrid(x);

z=-1*x.*y.*exp(-1*x.^2-y.^2); mesh(x,y,z)

6 作出椭球面x24?y2z2??1的图形. 91(该曲面的参数方程为

x?2sinucosv,y?3sinusinv,z?cosu, (0?u??,0?v?2?).)

>>

ezsurf('2*sin(u)*cos(v)','3*sin(u)*sin(v)','cos(u)',[0,pi,0,2*pi])

x2y2z27 作双叶双曲面2?2?2??1的图形.

1.51.41.3(曲面的参数方程是x?1.5cotucosv,y?1.4cotusinv,z?1.3cscu,其中参数

0?u??2,???v??时对应双叶双曲面的一叶, 参数??2?u?0,???v??时对应双叶

双曲面的另一叶.) >>

ezmesh('1.5*cot(u)*cos(v)','1.4*cot(u)*sin(v)','1.3*csc(u)', [0,pi/2,-pi,pi]); hold on

5

ezmesh('1.5*cot(u)*cos(v)','1.4*cot(u)*sin(v)','1.3*csc(u)',… [-pi/2,0,-1*pi,pi])

8 作出圆环

x?(8?3cosv)cosu,y?(8?3cosv)sinu,z?7sinv,(0?u?3?/2,?/2?v?2?)

的图形. >>

ezsurf('(8+3*cos(v))*cos(u)', …

'(8+3*cos(v))*sin(u)',…

'sin(v)',…

[0,3*pi/2,pi/2,2*pi])

9 作出球面x2?y2?z2?22和柱面(x?1)2?y2?1相交的图形. >>

[x,y,z]=sphere(50); surf(x,y,z) hold on

[x,y,z]=cylinder(1,50); x=x+1;

z(1,:)=-z(2,:); surf(x,y,z)

10 作出锥面x2?y2?z2和柱面(x?1)2?y2?1相交的图形. >>

ezsurf('u*cos(t)','u*sin(t)','u',[0,2*pi,-1,1]);

6

hold on

[x,y,z]=cylinder(1,50); x=x+1;

z(1,:)=-1; surf(x,y,z)

11用动画演示由曲线y?sinz,z?[0,?]绕z轴旋转产生旋转曲面的过程. （该曲线绕

z轴旋转所得旋转曲面的方程为x2?y2?sin2z, 其参数方程为

x?sinzcosu,y?sinzsinu,z?z,(z?[0,?],u?[0,2?])）

>>

m=moviein(1); for i=1:10

u=0:0.01:pi/5*i; z=0:0.01:pi;

[u,z]=meshgrid(u,z); x=sin(z).*cos(u); y=sin(z).*sin(u); z=z;

mesh(x,y,z);

m(:,i)=getframe; end

movie(m,1)

12.画出变上限函数?>>

syms xt

F=int(t*(sin(t))^2,0,x);

7

x0tsint2dt及其导函数的图形.

ezplot(F)

13.迪卡尔曲线

3at3at2x?,y?(x3?y3?3axy?0) 221?t1?t>>

ezplot('x^3+y^3-3*x*y')

at2at3x32,y?(y?) 14.蔓叶线x?a?x1?t21?t2>>

ezplot('y^2=x^3/(10-x)',[0,6,-6,6])

8

15.摆线x?a(t?sint),y?b(1?cost)。 >>

ezplot('y^2=x^3/(10-x)',[0,6,-6,6])

16.内摆线（星形线）x?acost,y?asint(x?y?a) >>

ezplot('(cos(t))^3','(sin(t))^3')

33232323

17.圆的渐伸线（渐开线)x?a(cost?tsint),y?a(sint?tcost) >>

ezplot('cos(t)+t*sin(t)','sin(t)-t*cos(t)')

18.空间螺线x?acost,y?bsint,z?ct >>

9

ezplot3('cos(t)','sin(t)','t')

19.阿基米德线r?a?,r?0。 >>

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

polar(r,t);

21.双纽线r2?a2cos2?((x2?y2)2?a2(x2?y2))。 >>

ezplot('(x^2+y^2)^2=30*(x^2-y^2)')

22.双纽线r2?a2sin2?((x2?y2)2?2a2xy)。 >>

ezplot('(x^2+y^2)^2=100*x*y')

10

(h)正螺面x?usinv,y?ucosv,z?4v >>

ezmesh('u*sin(v)','u*cos(v)','4*v');

16

实验二 一元函数微分学

1.分别画出坐标为(i,i2),(i2,4i2?i3),(i?1,2,?,10)的散点图, 并画出折线图. >>

i=1:10;

subplot(1,2,1) plot(i,i.^2,'r'); hold on;

plot(i,i.^2,'*'); subplot(1,2,2)

plot(i.^2,4*i.^2+i.^3,'r'); hold on;

plot(i.^2,4*i.^2+i.^3,'*');

2.画出前25个素数的散点图. >>

p=2:100;

for i=2:sqrt(100)

n=find(rem(p,i)==0&p~=i);

p(n)=[]; end

plot(p,'*

cosnx 4.讨论极限nlim??>> >> syms n %x=0时 syms n %x~=0时 f=(cos(0))^n; f=(cos(0))^n; limit(f,n,inf) limit(f,n,inf) >> >> >> >> ans = ans = 1 0. 所以limcosnx不存在。

n??

5. 在MATLAB中求下列极限（写出MATLAB命令和运行结果）

17

2sinx(1) lim(n?n?n) （2）lim(1?)3x （3） lim3

x??x?0x?3xn??x3x3?4x2?2(?1)n?4n2n3?1（4）lim （5）lim3. (6)limn?1 n?13n??5n?1x??n???3?47x?4sinx?xcosxlncotx(5) (7) lim (8) （9） limx2lnx lim2x??0x?0x??0lnxxsinx11?sinx? ?cos)x （12）lim??x??x??0x?0xx?x?>> 3/7 (10) limxx （11） lim(sin11?cosxsyms n x limit(sqrt(n+sqrt(n))-sqrt(n),inf) limit((1-2/x)^(3*x),inf) limit(sin(x)/(x^3+3*x),0) limit((3*x^3-4*x^2+2)/(7*x^3+4),inf) limit((2*n^3+1)/(5*n^3+1),inf) limit(((-1)^n+4^n)/(3^(n+1)+4^(n+1)),n,inf,'left') limit((sin(x)-x*cos(x))/(x^2*sin(x))) limit(log(cot(x))/log(x),x,0,'left') limit(x^2*log(x),x,0,'left') limit(x^x,x,0,'left') limit((sin(1/x)+cos(1/x))^x,inf) limit((sin(x)/x)^(1/(1-cos(x)))) >> >> ans =

1/2

ans = exp(-6)

ans = 1/3 ans =

6.讨论下列函数在指定点的连续性：

(1)函数?f(x)??x2?5x?6?,x??1在x??1处的连续性；?x???71x??1(2)函数?g(x)??1?x2?5x,x?0?在x??sinx?xx?00处的连续性； >>

18

ans = 2/5 ans = 1/4 ans = 1/3 ans = -1 ans = 0 ans = 1 ans = exp(1) ans = exp(-1/3)

syms x

f=(x^2-5*x-6)/(x+1); ff=-7;

if limit(f,x,-1)==ff

disp('f(x)?úx=-1′|á?D?'); end

gl=1-x^2-5*x; gr=sin(x)/x;

if limit(gl,x,0,'right')==limit(gr,x,0,'left') disp('g(x)?úx=0′|á?D?'); end >> >>

f(x)在x=-1处连续 g(x)在x=0处连续

7. 根据要求在MATLAB中求下列函数的导数

?1?x2?dyf(x)?arcsin?2???axaax1?xy?a?a?x?x??，求(1) ，求dx (2)

f??1???

d2y22??y?lnx?a?x22x?1y?xln(1?x)dydx（3）设，求 (4) ，求.

>> ans4=eval(ans4) syms x a >> >> y1=a^a+a^x+x^a+x^(a*x); ans1 = ans1=diff(y1) a^x*log(a)+x^a*a/x+x^(a*x)*(a*loy2=asin((1-x^2)/(1+x^2)); g(x)+a) ans2=diff(y2); ans2 = x='1'; -8.3264e-004 ans2=eval(ans2) ans3 = syms x [(1+1/(a^2+x^2)^(1/2)*x)/(x+(a^2y3=log(x+sqrt(a^2+x^2)); +x^2)^(1/2))dx] ans3=[diff(y3),'dx'] ans4 = y4=x^2*log(1+x); 10.7836 ans4=diff(y4,x,2); x='1';

??8.函数使

f(x)?1/x4在区间[1,2]上满足拉格朗日中值定理的条件, 因此存在??(1,2)f?(?)?f(f(2)?f(1))/(2?1).

可以验证这个结论的正确性. >> y=1/x.^4 x=2; y =

19

0.0625 x=1;

y=1/x.^4 y =1 syms x

diff(1/x.^4) ans = -4/x^5

x=1:0.01:2;

y=-4./x.^5;plot(x,y) 由如下图可知

9.验证拉格朗日定理对函数y?4x>> x=0;

y=4*x.^3-5*x.^2+x-2 y =-2 x=1;

y=4*x.^3-5*x.^2+x-2 y =-2

syms x;

diff(4*x.^3-5*x.^2+x-2) ans =

12*x^2 - 10*x + 1 x=0:0.01:1;

>> y=12*x.^2 - 10*x + 1;

plot(x,y)

由如右图可知成立：

10.证明:对函数y?px区间[a,b]的正中间. syms x;

syms p q r real; diff(p*x.^2+q*x+r)

20

23?5x2?x?2在区间[0,1]上的正确性.

?qx?r应用拉格朗日中值定理时, 所求得的点?总是位于

ans =

q + 2*p*x syms x;

syms p q r real;x=a; y=p*x.^2+q*x+r y =

p*a^2 + q*a + r syms x;

syms b p q r real;x=b; y=p*x.^2+q*x+r y =

p*b^2 + q*b + r syms x;

syms a b p q r real;x=(a+b)/2; z=q + 2*p*x z =

q + 2*p*(a/2 + b/2) syms x;

syms a b p q r real; y1=p*a^2 + q*a + r; y2=p*b^2 + q*b + r; z=(y2-y1)/(b-a) z =

(p*a^2 + q*a - p*b^2 - q*b)/(a - b) 经化简可知定理成立

d2y11.y?xln(1?x)，求2dx2x?1??

>> syms x

y=x^2*log(1+x); ans=diff(y,x,2); x='1';

ans=eval(ans) >> >> ans =

10.7836

?x?a(t?sint)dy12.设?，求

y?a(1?cost)dx?>>

21

syms t a

x=a*(t-sin(t)); y=a*(1-cos(t)); diff(y)/diff(x) >> >> ans =

sin(t)/(1-cos(t))

14x?2x3?3x?3613． 已知多项式f(x)?6x?2x?5x?1，，求：

（1）f(x)的根； (2) g(x)在闭区间[-1,2]上

532g(x)?的最小值；

（3）f(x)?g(x)，>>

f=[6,0,2,-5,0,1]; g=[0,1/6,2,0,-3,3]; ans1=roots(f) gg=polyder(g); x=roots(gg);

k=find(x<-1|x>2); x(k)=[];

ans2=min(polyval(g,x))

ans3_1=poly2sym(f+g) ans3_2=poly2sym(conv(f,g))

ans4=polyder(f) ans4=poly2sym(ans4) >> >> ans1 =

14. 已知函数

f(x)f(x)?g(x)和g(x)； （4）f(x)的导数。

-0.4056 + 0.9685i -0.4056 - 0.9685i 0.6065 + 0.0915i 0.6065 - 0.0915i -0.4018 ans2 = 1.6245 ans3_1 =

6*x^5+1/6*x^4+4*x^3-5*x^2-3*x+4 ans3_2 =

x^9+12*x^8+1/3*x^7-89/6*x^6+8*x^5-35/6*x^4+23*x^3-15*x^2-3*x+3 ans4 =

30 0 6 -10 0 ans4 =

30*x^4+6*x^2-10*x

f(x)?16254x?2x5?x?60x3?150x2?180x?25, 22在区间[?6,6]上画出函数f(x),f?(x),f??(x)的图形, 并找出所有的驻点和拐点.

22

>>

syms x real

f=1/2*x^6-2*x^5-25/4*x^4+60*x^3-150*x^2-180*x-25;

ff=diff(f); fff=diff(f,2); subplot(2,1,1) ezplot(f,[-6,6]); subplot(2,2,3) ezplot(ff,[-6,6]); subplot(2,2,4)

ezplot(fff,[-6,6]);

ans1=eval(solve(ff)); %求驻点

k=find(imag(ans1)~=0|ans1<-6|ans1>6); ans1(k)=[]

ans2=eval(solve(fff)); %求拐点

k=find(imag(ans2)~=0|ans2<-6|ans2>6); ans2(k)=[] >> >> ans1 =

3.5367 -0.4642 -3.9608 ans2 =

1.3417 -2.8623

23

x?15.求函数y?2sin2(2x)?5xcos2???的位于区间(0,?)内的极值的近似值.

2?2?先建立函数文件f.m : function y=f(x)

y=2*sin(2*x)*sin(2*x)+5/2*x*cos(x/2)*cos(x/2); >>

x=fminbnd('-f(x)',0,1) %求极大值点 f(x) %求极大值 >> >> x =

0.8642 ans =

3.7323 >>

x=fminbnd('f(x)',0.5,2.5) %求极小值点 f(x) %求极小值点 >> >> x =

1.6239 ans =

1.9446 >>

x=fminbnd('-f(x)',1.5,pi) %求极大值点 f(x) %求极大值点 >> >> x =

2.2449 ans =

2.9571

实验三 一元函数积分学

1．用MATLAB计算下列不定积分。 （1）?x2?1dx （2）?axsinxcos2xdx 2x>>

syms x a

ans1=int(sqrt(x^2+1)/x^2)

ans2=int(a^x*sin(x)*(cos(x))^2,x) 2．用MATLAB求解下列各积分。 （1）?

2?0?ecosxdx （2）?e?tsin2tdt

02x24

2?x20?x?1 （3）设f(x)??，求?f(x)dx。

0?x1?x?2>>

syms x t

ans1=int(exp(2*x)*cos(x),0,2) ans2=int(exp(-t)*sin(2*t),0,inf) ans3=int(x^2,0,1)+int(x,1,2) >> >>

ans1 =2/5*exp(4)*cos(2)+1/5*exp(4)*sin(2)-2/5 ans2 =2/5 ans3 =11/6

4．求由曲线x2?(y?5)2?16绕x轴旋转所产生的旋转体的体积。

>> syms x

V=eval(int(pi*(sqrt(16-x^2)+5)^2,-4,4)-int(pi*(sqrt(16-x^2)-5)^2,-4,4)) >> >>

V =1.5791e+003

5．求下列曲线与所围成图形的面积： （1）y?12； （2）r?2sin?与r2?cos2? x与x2?y2?8（两部分都要计算）

2>>

syms x y %求(1)

S1=int(1/2*x^2+sqrt(8-x^2),-2,2) S2=int(sqrt(8-x^2)-1/2*x^2,-2,2) >> >> S1 =12.9499 S2 =7.6165

>>

syms t %求（2） r1=sqrt(2)*sin(t); r2=sqrt(cos(2*t));

S=2*double(int(int(1,r1,r2),0,pi/6)) >> >> S =0.5691

25

x26.计算半立方抛物线y2?(x?1)3被抛物线y2?截得的一段弧的长度。

33>>% 注意利用对称性 syms x

y=sqrt(2/3*(x-1)^3);

l=2*double(int(sqrt(1+(diff(y))^2),1,2)) >> >> l =2.6248

26

实验四 多元函数微积分

求多元函数的偏导数与全微分

1.1>>

syms x y;z=sin(x*y)+cos(x*y)^2; zx=diff(z,x) zy=diff(z,y) zzxx=diff(z,x,2) zzxy=diff(zx,y) >> >>

zx =cos(x*y)*y-2*cos(x*y)*sin(x*y)*y zy =cos(x*y)*x-2*cos(x*y)*sin(x*y)*x

zzxx =-sin(x*y)*y^2+2*sin(x*y)^2*y^2-2*cos(x*y)^2*y^2

zzxy=-sin(x*y)*x*y+cos(x*y)+2*sin(x*y)^2*x*y-2*cos(x*y)^2*x*y-2*cos(x*y)*sin(x*y)

1.2>>

syms x y u v; vx=-fx/f1v f1=exp(u)+u*sin(v)-x; vy=-fy/f2v f2=exp(u)-u*cos(v)-y; >> >> f1u=diff(f1,u); ux = f1v=diff(f1,v); 1/(exp(u)+sin(v)) fx=diff(f1,x); uy = f2u=diff(f2,u); 1/(exp(u)-cos(v)) f2v=diff(f2,v); vx = fy=diff(f2,y); 1/u/cos(v) ux=-fx/f1u vy = uy=-fy/f2u 1/u/sin(v) 微分学的几何应用

1.3 求出曲面z?2x2?y2在点(1,1)处的切平面、法线方程, 并画出图形. >>

[x,y]=meshgrid(-5:0.1:5); z=2.*x.^2+y.^2; mesh(x,y,z) hold on

[x,y]=meshgrid(-10:0.1:10);z=4*x+2*y-3; plot3(x,y,z) hold on

line([41,-39],[21,-19],[-7,13]) axis([-20 20 -20 20 -40 40])

27

?z?z?2z?2z,,2,.2z?sin(xy)?cos(xy),?x?y?x?y?x设求

?u?u?v?v,,,.x?e?usinv,y?e?ucosv?x?y?x?y 设,求

uu

1.4求曲面k(x,y)?4?1164?在点?,,?x2?y2?1?4221?处的切平面方程, 并把曲面和它的切平

面作在同一图形里.

>>

syms x y k;

df_dx=diff(4/(x^2+y^2+1),x) df_dy=diff(4/(x^2+y^2+1),y) a=linspace(-10,10,100); b=a;

[a,b]=meshgrid(a,b); c=4./(a.^2+b.^2+1);

d=-8/((1/4)^2+(1/2)^2+1)^2*(1/4);

e=-8/((1/4)^2+(1/2)^2+1)^2*(1/2);

f=d.*(a-1/4)+e.*(b-1/2)+64/21; mesh(a,b,c); hold on;

mesh(a,b,f);

axis([-10,10,-10,10,-2,5]);

多元函数的极值

1.5求f(x,y)?x3?y3?3x2?3y2?9x的极值. >>

syms x y

f=x^3-y^3+3*x^2+3*y^2-9*x; fx=diff(f,x); fy=diff(f,y);

[x y]=solve(fx,fy) %求极值点

>> >>

28

x=1 -3 y=0 2 >>

x=1;y=0; eval(f) x=-3;y=2; eval(f)

>> >> ans = -5 ans = 31

1.6 求函数z?x2?y2在条件x2?y2?x?y?1?0下的极值. >>

syms x y

y=solve('x^2+y^2+x+y-1',y) %把y用x表示出来 >> >> y =

-1/2+1/2*(5-4*x^2-4*x)^(1/2) -1/2-1/2*(5-4*x^2-4*x)^(1/2) >>

z=x^2+(-1/2+1/2*(5-4*x^2-4*x)^(1/2))^2; %将z用x表示出来 ezplot(z) zx=diff(z,x);

x=eval(solve(diff(z))); %求出极值点 eval(z) %求极值 >> >> ans =

0.2679 >>

syms x y

z=x^2+(-1/2-1/2*(5-4*x^2-4*x)^(1/2))^2; %将z用x表示出来 ezplot(z) zx=diff(z,x);

x=eval(solve(diff(z))); &求极值点 eval(z) %求极值 >> >> ans =

3.7321

实验2 多元函数积分学（基础实验） 计算重积分

2.1计算??xy2dxdy, 其中D为由x?y?2,x?Dy, y?2所围成的有界区域.

29

>>

syms x y;

int(int(x*y^2,x,2-y,y^0.5),y,1,2) >> >> ans =

193/120

2.2计算???(x2?y2?z)dxdydz, 其中?由曲面z??2?x2?y2与z?x2?y2围成.

>> syms t r z;int(int(int((r^2+z)*r,z,r,(2-r^2)^0.5),r,0,1),t,0,2*pi) >> >> ans =

2.1211

重积分的应用

2.3 求由曲面f?x,y??1?x?y与g?x,y??2?x2?y2所围成的空间区域?的体积. >>

syms t r;

int(int((3/2-r^2)*r,r,0,(3/2)^0.5),t,0,2*pi) >> >> ans =

3.5343

2.4 在Oxz平面内有一个半径为2的圆, 它与z轴在原点O相切, 求它绕z轴旋转一周所得旋转体体积. >>

syms x;

int(4*pi*x*(4-(x-2)^2)^0.5,x,0,4) >> >> ans =

157.9137

计算曲线积分

2.5求 ?f(x,y,z)ds, 其中f?x,y,z??1?30x2?10y,积分路径为

L L:x?t,y?t2,z?3t2,0?y?2.

(注意到,弧长微元ds?xt2?yt2?zt2dt, 将曲线积分化为定积分) >>

syms t; x=t; y=t^2; z=3*t^2;

f=diff([x,y,z],t);

30

fun=inline('((1+30*t.^2).^0.5+10*t.^2).*(1+40*t.^2).^0.5','t'); quad(fun,0,2) >> >> ans =

348.9428

2.6求?LF.dr, 其中

F?xy6i?3x(xy5?2)j,r(t)?2costi?sintj,0?t?2?.

>>

syms t; x=cos(t); y=sin(t);

int(x*y^6*(-2*sin(t))+3*x*(x*y^5+2)*cos(t),t,0,2*pi) >> >> ans = 6*pi

计算曲面积分

2.7计算曲面积分??(xy?yz?zx)dS, 其中?为锥面z?x2?y2被柱面x2?y2?2x所截

?得的有限部分.

(注意到,面积微元dS?21?zx?z2ydxdy, 投影曲线x2?y2?2x的极坐标方程为

?2?t?r?2cost,??2,

将曲面积分化作二重积分,并采用极坐标计算重积分.)

>>

syms t r; x=r*cos(t); y=r*sin(t); z=r;

int(int((x*y+y*z+x*z)*r*2^0.5,r,0,2*cos(t)),t,-pi/2,pi/2) >> >> ans =

6.0340

2.8计算曲面积分??x3dydz?y3dzdx?z3dxdy, 其中?为球面x2?y2?x2?a2的外侧.

?>>syms t s r; syms a real;

int(int(int(3*r^4*sin(s),r,0,a),s,0,pi),t,0,2*pi) >> >> ans =

12/5*a^5*pi

31

曲线拟合

3.1 为研究某一化学反应过程中温度x(?C)对产品得率y(%)的影响, 测得数据如下:

x 100 110 120 130 140 150 160 170 180 190 y 45 51 54 61 66 70 74 78 85 89

试求其拟合曲线. >>

x=[100,110,120,130,140,150,160,170,180,190]; y=[45,51,54,61,66,70,74,78,85,89]; a=polyfit(x,y,1) >> >>

a = 0.4830 -2.7394

即拟合曲线为：y=0.4830x-2.7394

3.2 给定平面上点的坐标如下表:

x0.10.20.30.40.50.60.70.80.9

y5.12345.30575.56875.93786.43377.09787.94939.025310.3627试求其拟合曲线.

>>

x=[0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9];

y=[5.1234,5.3057,5.5687,5.9378,6.4337,7.0978,7.9493,9.0253,10.3627]; a=polyfit(x,y,3) >> >>

a = 4.9875 0.6902 1.3202 4.9774

即拟合曲线为：y=4.9875x^3+0.6902x^2+1.3202x+4.9774

3.3已知lg(x)的[1，101]区间11个整数采样点的函数值如下表

511.7076 试求lg(x)的5次拟合多相式p(x)，并分别绘制出lg(x)和p(x)在[1，101]区间的函数曲线 >>

x=[100,110,120,130,140,150,160,170,180,190]; y=[45,51,54,61,66,70,74,78,85,89]; p=polyfit(x,y,5) >> >> p =

1.0e+014 *

-0.0000 0.0004 -0.0226 0.5192 -4.5093 4.0122

x lg(x)10111.0414211.3222311.4914411.6128

32

实验六 无穷级数及微分方程 （基础实验）

数项级数

1.1(1) 观察级数??n?11的部分和序列的变化趋势. 2n?1nn? (2)观察级数?1的部分和序列的变化趋势. (1) x=0;

for n=1:50; x=x+1/n^2;

plot(n,x,’r*’) hold on end

(2) x=0;

for n=1:100; x=x+1/n;

plot(n,x,'r*') hold on end

33

10n1.2 设an?, 求

n!?an?1?n.

>> s=10;

for i=1:inf;

s=s+s*10/(i+1); if s*10/(i+1)==eps break; end end s

>> >> s =

5.2257e+086

求幂级数的收敛域 1.3 求?n?0?42n(x?3)nn?1的收敛域与和函数.

>>

syms n x

limit(((4^(2*(n+1))*(x-3)^(n+1))/(n+1+1))/((4^(2*n)*(x-3)^n)/(n+1)),n,inf) >> >> ans = 16*x-48 >>

syms n x

x=solve('16*x-48=1') x=solve('16*x-48=-1') >> >> x = 49/16 x = 47/16

所以收敛域为（47/16,49/16) >>

syms n x

symsum(4^(2*n)*(x-3)^n/(n+1),n,1,inf) >> >> ans =

-1/16/(x-3)*(16*x-48+log(49-16*x))

?3n 1.4 设an?, 求?ann!n?1.

34

>> s=3;

for i=1:inf;

s=s+s*3/(i+1); if s*3/(i+1)==eps break; end end s

>> >> s =

1.2379e+027

1.5求下列级数的和：

I1??n?1?2n?12n，I2??n?1?12(2n?1)，I3??n?1?(?1)n?1n，I4??n2

n?1100>> syms n

I1=symsum((2*n-1)/2^n,n,1,inf) I2=symsum(1/2/(2*n+1),n,1,inf) I3=symsum((-1)^(n+1)/n,n,1,inf) I4=symsum(n^2,n,1,100)

求幂级数的收敛域 1.6 求?(?1)n?1?n?1>> >> I1 =3 I2 =Inf I3 =log(2) I4 =338350

xnn的收敛域与和函数.

>>

syms n x

limit(((-1)^(n-1+1)*x^(n+1)/(n+1))/((-1)^(n-1)*x^n/n),n,inf) >> >> ans = -x

所以收敛域为（-1,1） >>

syms n x

limit(((-1)^(n-1)*(x)^n/n),n,inf) >> >> syms n x

limit((-1)^(n-1)*(x)^n/n,n,inf) >> >> ans =

limit((-1)^(n-1)*x^n/n,n = Inf)

函数的幂级数展开

1.7将函数sinx展开为x的幂级数，分别展开至5次和20次。

35

>> syms x

ans1=taylor(sin(x),x,6,0) ans2=taylor(sin(x),x,21,0) >> >> ans1 =

x-1/6*x^3+1/120*x^5 ans2 =

x-1/6*x^3+1/120*x^5-1/5040*x^7+1/362880*x^9-1/39916800*x^11+1/6227020800*x^13-1/1307674368000*x^15+1/355687428096000*x^17-1/121645100408832000*x^19

1.8将函数(1?x)m展开为x的幂级数，m为任意常数。展开至4次幂。 >>

syms x m

taylor((1+x)^m,x,5,0) >> >> ans =

1+m*x+1/2*m*(m-1)*x^2+1/6*m*(m-1)*(m-2)*x^3+1/24*m*(m-1)*(m-2)*(m-3)*x^4

11.9将函数f(x)?2展开为(x?2)的幂级数。

x?5x?3>>

syms x m

taylor(1/(x^2+5*x-3),x,5,2) >> >> ans =

29/121-9/121*x+70/1331*(x-2)^2-531/14641*(x-2)^3+4009/161051*(x-2)^4

?1.10将函数cosx展开成(x?)的幂级数，取前10项。

3>> syms x

taylor(cos(x),x,10,pi/3) >> >> ans =

1/2-1/2*3^(1/2)*(x-1/3*pi)-1/4*(x-1/3*pi)^2+1/12*3^(1/2)*(x-1/3*pi)^3+1/48*(x-1/3*pi)^4-1/240*3^(1/2)*(x-1/3*pi)^5-1/1440*(x-1/3*pi)^6+1/10080*3^(1/2)*(x-1/3*pi)^7+1/80640*(x-1/3*pi)^8-1/725760*3^(1/2)*(x-1/3*pi)^9

1.11求函数f(x)?x2在[??,?]上的傅立叶级数。 >> syms x f=x^2

a0=fa(f,0)

36

n=1:5; a=fa(f,n) b=fb(f,n) >> >> a0 =

2/3*pi^2 a =

[ -4, 1, -4/9, 1/4, -4/25] b =

[ 0, 0, 0, 0, 0] 所以fx=2?2 414?4cosx?1cos2x?cos3x?cos4x?cos5x?...394251.12求出函数f(x)?x3?x2在区间[??,?]上的前11个傅立叶系数，即n＝5。 >>

function an=fa(f,n) %求an syms x

an=int(f*cos(n*x),x,-pi,pi)/pi; function bn=fb(f,n) % 求bn syms x

bn=int(f*sin(n*x),x,-pi,pi)/pi; >> >> a0 =

2/3*pi^2 n =

1 2 3 4 5 a =

[ -4, 1, -4/9, 1/4, -4/25] b =

[ (2*pi^3-12*pi)/pi, (-pi^3+3/2*pi)/pi, (2/3*pi^3-4/9*pi)/pi, (-1/2*pi^3+3/16*pi)/pi, (2/5*pi^3-12/125*pi)/pi]

1.13求arctanx的5阶泰勒展开式. >>

syms x;taylor(atan(x)) >> >> ans =

x-1/3*x^3+1/5*x^5

221.14 求e??x?1??x?1?在x?1处的8阶泰勒展开, 并通过作图比较函数和它的近似多 >>

syms x;

t=taylor(exp(-(x-1)^2*(x+1)^2),9,1)

37

ezplot(t,[0,2]); hold on;

ezplot(exp(-(x-1)^2*(x+1)^2),[0,2]) >> >> ans =

1-4*(x-1)^2-4*(x-1)^3+7*(x-1)^4+16*(x-1)^5+4/3*(x-1)^6-28*(x-1)^7-173/6*(x-1)^8

1.15 利用泰勒公式ex?1?x?x|Rn|?0.005,问

22!???xn!n

?Rn(x)近似计算ex. 若|x|?1,要求截断误差

n应取多大?

1.16 观察函数f(x)?sinx各阶泰勒展开的图形. (1) 固定x0?0,观察阶数n的影响;

>> syms x

f10=ezplot(taylor(sin(x),x,10,0)); set(f10,'color','r') hold on

f50=ezplot(taylor(sin(x),x,50,0)); set(f50,'color','b') hold on

f100=ezplot(taylor(sin(x),x,100,0)); set(f100,'color','c') title('sin x')

(2) 扩大显示区间范围, 以观察在偏离展开点x0时泰勒多项式对函数的逼近情况;

(3) 固定n?10,观察x0的影响.

38

syms x

f50=ezplot(taylor(sin(x),x,50,1)); set(f50,'color','r') hold on

f50=ezplot(taylor(sin(x),x,50,5)); set(f50,'color','b') hold on

f50=ezplot(taylor(sin(x),x,50,10)); set(f50,'color','c') title('sin x')

求解微分方程

2.1求微分方程 y??2xy?xe?x的通解. >>

y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') >> >>

y =(1/2*x^2+C1)*exp(-x^2)

2.2求微分方程xy??y?ex?0在初始条件yx?1?2e下的特解.

2>>

y=dsolve('x*Dy+y=exp(x)','y(1)=2*exp(1)','x')

>> >>

y =(exp(x)+exp(1))/x

2.3求解微分方程y???2x?ex, 并作出其积分曲线. >>

y=dsolve('D2y-2*x=exp(x)','x') >> >> y =

1/3*x^3+exp(x)+C1*x+C2 >>

ezplot('1/3*x^3+exp(x)+1*x+2')

39

2.4

?dxt?dt?x?2y?e求微分方程组?dy??x?y?0??dt在初始条件xt?0?1,yt?0?0下的特解.

>>

[x y]=dsolve('Dx+x+2*y-exp(t)','Dy-x-y','x(0)=1','y(0)=0','t') >> >>

x =cos(t)

y =1/2*sin(t)-1/2*cos(t)+1/2*exp(t)

2.5求出初值问题

22??y???y?sinx?y?cosx ???y(0)?1,y(0)?0?的数值解, 并作出数值解的图形.

>>

function dy=ffer(x,y) dy=zeros(2,1);dy(1)=y(2);

dy(2)=-y(2)*(sin(x))^2-y(1)+(cos(x))^2; [X,Y]=ode23s('ffer',[0 4],[1 0]) plot(X,Y(:,1),'-') >> >>X = 0.9767 0 1.1374 0.0183 1.2859 0.0366 1.4337 0.0550 1.5851 0.0733 1.7415 0.0917 1.9018 0.1100 2.0632 0.1284 2.2229 0.1469 2.3794 0.1657 2.5324 0.1870 2.6836 0.2111 2.8377 0.2384 3.0051 0.2694 3.1122 0.3048 3.2193 0.3453 3.3242 0.3920 3.4361 0.4461 3.5540 0.5098 3.6787 0.5863 3.8127 0.6823 3.9528 0.8161 4.0000

40

本文来源：https://www.bwwdw.com/article/9bd7.html