实验六常微分方程

函数表达Matlab表达函数表达Matlab表达

y=sinxy=sin(x)y=exy=exp(x)

y=cosxy=cos(x)y=Inxy=iog(x)

y=x2y=xA2y=arctanxy=atan(x)

sinx

y=x2sinxy=x2*sin(x)y=xy=sin(x)/x

实验六常微分方程

dsolve]dx=-a*x')%输出ans=Cl*exp(-(2*x)

1=/%{'小=一〃*工；比(0)=1'，'5)%输出x=exp(-〃*s)

w=dsolve(D3w=-w7w(0)=1,Dw(0)=0,D2w(0)=0,)%输出略

[7,g]=^Zve('zy=/+g,Dg=-/+g','/(o)=r,'g(o)=2')%输出略

一、实验内容

求微分方程的解析解

2,

例1求微分方程y+2.二比一"的通解。

输入：

dsolve(*dy+2*x*y=x*exp(-xA2)',*x*)

便得到微分方程的通解：

(l/2*xA2+Cl)*exp(-x-'2)

其中Cl是任意常数。

例2求微分方程=0在初始条件y(l)=2e下的特解。

输入：

dsolve(*x*dy+y-exp(-x)=0*z*y(1)=2*exp(1)*z*x*)

便得到微分方程的特解：

(-exp(-x)+exp(-1)+2*exp(1))/x

实验七空间图形的画法

(一)、三维曲线的绘制

(1)plot3命令

绘制三维曲线，基本使用形式是：

plot3(x,y,z,'s')

例如，一条空间螺旋线的参数方程是：x=cosf,y=sinf,z=看(0W8%)。

输入：

t=0:0.1:8*pi;

x=cos(t);

y=sin(t);

z=t/10;

plot3(x,y,z)

xlabel(1x1);

ylabel(1y1);

zlabel(1z1);

则输出了一条螺旋线(见图6-1)o

图6-1

(2)ezplot3命令

同绘制二维曲线相似，也有简洁的绘制空间曲线命令ez〃"3,具体使用方法

同相似，形式是：

ezpI3o(fx(t)*,*y(t)*,*z

其中x(t),y(t),z(t)是曲线的参数方程的表达式。tl,t2是作图时参数t的范围。

例如上述螺线也可输入下面命令得到：

/p，“3('cos(t)；'sin(t)7t/10\[0,8*pi\)

输出图形略。

(二)三维曲面网线图与曲面图的绘制

Meshgrid：此函数用来产生三维绘图时的阵列(网格).

例：

x=l:3;

y=10:14;

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

X=

123

123

123

123

123

101010

111111

121212

131313

141414

mesh:绘制3d网状立体图

surf:绘制3d彩色表面图

title:标识文字，title('text',…)

xlabel：对x轴的坐标

ylabel：对y轴的坐标

gridon:在图表中交网格线

holdon:保持目前的图表，后续可附加

subplot(m,n,p)：多个窗口同时显示，子窗口图形个数为m*n矩阵

MATLAB软件绘制曲面图基本形式是：

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

(2)Z=/(X,Y)

(3)mesh(X,Y,Z)%绘制网线图

(4)5wrf(X,Y,Z)%绘制曲面图

例如，画出曲面Z=%2+y2的图形。输入：

x=-2:0.1:2;

y=-2:0.1:2;

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

Z=X.A2+Y.A2;

Mesh(X,Y,Z)

surf(X,Y,Z)

得到曲面Z=x2+y2(见图6-2)。

图6-2

又如，参数方程：x=2sin℃ose,y=2sin°sine,z=2cos0是以原点为圆心,

半径为2的球面，其中0<0<万万，因此只要输入：

t=0:0.1:pi;

r=0:0.1:2*pi;

[R,T]=meshgrid(r,t);

x=2.*sin(T).*cos(R);

y=2.*sin(T).*sin(R);

z=2,*cos(T);

surf(x,y,z)

便作出了方程为4+y2+z2=4的球面(见图6・3方

图6-3

绘制二元函数图形也可用简捷绘制的ezs〃^命令，它的使用格式为:

ezsurf(f(x,y),[a,b,u,v])

即可绘制函数/(x,y)在区域［a,勿上的图形。当省略区域时，默认区域

是[-2肛21]x[-2乃,21]o

一、实验内容

(一)空间曲线

例1作出空间曲线x=/cos/,y=/sin/,z=2/(0W/W6万)的图形。

输入：

ezplot3('t*cos(t)'；t*sin(t)','2*t',[0,6*pi])

输出如图6-4所示。输出如图6-4所示。

（二）空间二次曲面

2,2

yz，

例2作椭球面工+—+—=1的图形。

491

该曲面的参数方程是x=2sin^>cos0,y=3sin^sin0,z=cos(p,其中

G<(p<7i.O<0<o输入：

t=0:0.1:pi;

r=0:0.1:2*pi;

[R,T]=meshgrid(r,t);

x=2,*sin(T),*cos(R);

y=3.*sin(T).*sin(R);

z=cos(T);

surf(x,y,z)

输出如图6-5所示。

图6-5

222

例3作单叶双曲面三+匕—二=1的图形。

149

该曲面的参数方程是％=se&・svny=,2«secvo,其中

jrjr

——<w<—,Q^v<TQ.O输入：

22

t=-pi/4:0.1:pi/4;

r=0:0.1:2*pi;

[R,T]=meshgrid(r,t);

x=sin(R).*sec(T);

y=cos(R).*sec(T);

z=3*tan(T);

surf(x,y,z)

输出如图6-6所示。

图6-6

222

例4作双叶双曲面三+三一三=一1的图形。

该曲面的参数方程是x=1.5cotw-cosv,=1.4cotw-sinv,z=1.3cscw,其中参数

rr.TT

0<M«—乃4vW万对应双叶双曲面的一叶，参数——〈”<0,—乃乃对应双

22

叶双曲面的另一叶。输入：

t=pi/1000:0.1:pi/2;

r=-pi:0.1:pi;

[R,T]=meshgrid(r,t);

xl=l.5.*cos(R),*cot(T);

yl=l.4.*sin(R).*cot(T);

zl=l.3.*csc(T);

mesh(xl,yl,zl)

holdon

mesh(-xlr-yl,-zl)

输出如图6-7所示。

例3求微分方程产2)/+5y=e，cos2x的通解。(simplify)化简结果

输入：

simplify(dsolve(*D2y-2dy+5*y=exp(x)*cos(2*x)*x*))

便得到微分方程的通解：

sin(5A(1/2)*x)*C2+cos(5A(1/2)*x)*Cl+2/5+l/10*exp(x)*cos(2*x)

+1/5*exp(x)*sin(2*x)

其中Cl,C2是任意常数。

解微分方程组的命令格式为：

[x,y]=dsolve^dx=/(x,y)'6=g(x,y)1)

dxct

---bx+2y=e黑))::下的特解。

例4求微分方程组出在初始条件＜

电-x-y=0

〔dt

输入:

[xzy]=dsolve('dx=-x-2*y+exp(t)','dy=x+y','x(0)=l','y(0)=0')

便得到微分方程的特解:

x=cos(t)

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

实验八

(-)求偏导数命令

命令d疗既可以用于求一元函数的导数，也可以用于求多元函数的偏导数。

用于求偏导数时，可根据需要分别采用如下几种形式：

y,z),x)%求/(%,y/)对刀的偏导数

或「(/(X,Xz),y)%求于(x,y,z)对y的偏导数

d疗(/(x,y,z),无,2)或diff(diff(f(x,y,z),x),x)%求/(x,y,z)对％的

二阶偏导数

diff{diff(f(x,y,z),x),y%求/(%y,z)对x、y的混合偏导数

其余类推。

一、实验内容

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

例1iSz=sin(xy)+cos2(xy),求当

oxoyoxoxoy

输入：

symsxy

z=*sin(x*y+(cos(x*y))A2)*;

diff(z,x)

diff(z,y)

diff(z,x,2)

diff(diff(z,x),y)

便依次得到函数表达式及所求的四个偏导数的结果：

ans=cos(cos(x*y)A2+x*y)*(y-2*y*cos(x*y)*sin(x*y))

ans=cos(cos(x*y)A2+x*y)*(x-2*x*cos(x*y)*sin(x*y))

ans=-sin(cos(x*y)A2+x*y)*(y-2*y*cos(x*y)*sin(x*y))A2-

cos(cos(x*y)A2+x*y)*(2*yA2*cos(x*y)A2-2*yA2*sin(x*y)A2)

ans=-cos(cos(x*y)A2+x*y)*(2*x*y*cos(x*y)A2+

2*cos(x*y)*sin(x*y)-2*x*y*sin(x*y)A2-1)-sin(cos(x*y)A2+x*y)*(x

-2*x*cos(x*y)*sin(x*y))*(y-2*y*cos(x*y)*sin(x*y))

例2设z=(a+砂)"求二，一。

dxdy

输入：

symsxy

z=*(a+x*y)Ay*;

diff(z,x)

diff(z,y)

则有输出：

ans=yA2*(a+x*y)A(y-1)

ans=log(a+x*y)*(a+x*y)Ay+x*y*(a+x*y)A(y-1)

(二)多元函数的极值

例3求函数/(%»)=%3-/+3_¥2+3,2一9%的极值。

输入：

symsxy

f=,xA3-yA3+3*xA2+3*yA2-9*x,;

fx=diff(f,x)

fy=diff(f,y)

输出为：

fx=3*xA2+6*x-9

fy=6*y-3*yA2

再输入：

[x0,y0]=solve('3*x/'2+6*x-9=0','-3*y/'2+6*y','x','y')

输出为如下四个驻点：

xO=

1

-3

1

-3

yO=

0

0

2

2

再输入：

fxx=diff(f,x,2);

fyy=diff(f,y,2);

fxy=diff(fx,y);

a=(fxx)*(fyy)-(fxy)A2;

x=-3;

y=0；

al=eval(a)

bl=eval(fxx)

cl=eval(f)

x=-3;

y=2;

a2=eval(a)

b2=eval(fxx)

c2=eval(f)

x=l;

y=0;

a3=eval(a)

b3=eval(fxx)

c3=eval(f)

x=l;

y=2;

a4=eval(a)

b4=eval(fxx)

c4=eval(f)

我们得到了四个驻点处的判别式函数A(x,y)=%4-抬，儿与/的值。归

纳以后用表格形式列出。

y九A(x,y)f

-30-12-7227

-32-127231

101272-5

1212-72-1

从中可见：

%=_3,y=2时，判别式A(x,y)=72,九=-12,因此函数有极大值31.

x=l,y=O时，判别式4(x,y)=72,九=12,因此函数有极小值-5.

龙=—3,y=0和无=l,y=2时，判别式A(x,y)=—72,函数在这些点不取极值.

二、实验作业

“、儿2_p.&dz

1.设z=e”,求——,—o

dxdy

-(x2+j;2)//822

2.z=^(cosx+siny),求包2,。％o

dxdydxdy

3.求f(x,y)=-120x3-30x4+18x5+5x6+30xy2的极值。

实验八

实验目的

掌握用MATLAB计算二重积分的方法。提高应用重积分解决各种应用问题的

能力。

学习MATLAB命令

(-)重积分命令

命令int也可用于计算重积分，常用格式如下：

int(int(/(x,y),y,%(x),y2(x),x,a,b)

用来求二重积分力力y,其中。是X型区域aWxWA,%(x)WyW%(x)。

D

注：y型区域上的二重积分的格式也是类似的。

例如要计算xy2dydx,

输入：

symsxy

int(int(x*yA2,0,x),x,O,l)

得到：

ans=1/15

二、实验内容

(一)计算重积分

例1计算口划其中。为由x+y=2,x=6，y=2所围成的有界

D

区域。

先做出区域。的草图，手工就可以确定积分限。应先对X积分，

输入：

symsxy

A

int(int(x*y2zx,2-y,sqrt(y)),y,l,2)

输出为：

ans=193/120

例2计算口"(,+廿)办dy,其中。为Y+V41。

D

如果用直角坐标计算，输入：

symsxyreal

f=exp(-(xA2+yA2));

int(int(f,y,-sqrt(l-xA2),sqrt(l-xA2)),x,-l,l)

输出为:

ans=int((piA(1/2)*erf((1-xA2)A(1/2)))/exp(xA2),x=-1..1)

积分遇到了困难。改用极坐标，也用手工确定积分限，输入:

symsrs

f=exp(-(r人2))*r;

int(int(f,r,0,1),s,0,2*pi)

输出为:

ans=-pi*(1/exp(1)-1)

（二）重积分的应用

例3求曲面z=4-x?-y2在My平面上的面积s。

输入：

symsxyz

ezsurf(14-xA2-yA21)

输出如图8-1所示。

图8-1

观察曲面的图形，可见是一个旋转抛物面。计算曲面面积公式是

S=JJyjl+z[+z^dxdy。

D

输入：

symsxy

z=*4-xA2-yA2,;

f=sqrt(1+diff(z,x)A2+diff(z,y)人2)

输出为：

f=(4*x-2+4*yb+1)A(1/2)

因此用极坐标计算。

输入：

symsrt

f=*sqrt(l+4*rA2)*r*;

s=int(int(f,r,0,2),t,0,2*pi)

输出为：

s=(pi*(17*17A(1/2)-1))/6

三、实验作业

1.计算［%jj（ysin九一xsin丁）为协。

2.求积分Jsin（exy）dxdy的近似值。

用极坐标计算££2：2dydxo

实验九二重积分

（二）重积分命令

命令int也可用于计算重积分，常用格式如下:

int(int(/(x,y),y,%(x),y2(x),x,a,b)

用来求二重积分U/(%,y)dxdy,其中。是X型区域。4%<九%(%)<%(》)。

D

注：F型区域上的二重积分的格式也是类似的。

例如要计算xy2dydx,

symsxy

int(int(x*yA2,y,0,x),x,O,l)

得到：ans=1/15

(三)二元函数的数值积分

函数加的功能是求矩形区域上的二元函数的数值积分。其格式如下:

q=dblquad{jun,xmin,xmax,ymin,ymax,tol)

最后一个参数的意义是用指定的精度⑹代替默认精度KF',再进行计算。

例如输入：

f=inline(*sqrt(max(1-(x.A2+y.A2),0)),);

Q=dblquad(f,-1,1,-lf1)

得到：

Q=

2.0944

四、实验内容

计算重积分

例1计算孙2dmy,其中。为由x+y=2,x=6，y=2所围成的有界区域。

D

先做出区域。的草图，手工就可以确定积分限。应先对X积分，输入：

symsxy

A

int(int(x*y2zx,2-y,sqrt(y)),y,l,2)

输出为：

ans=193/120

例2计算JJ如力，其中。为好+J?<i。

D

如果用直角坐标计算，输入:

symsxyreal

f=exp(-(xA2+yA2));

int(int(f,y,-sqrt(l-xA2),sqrt(l-xA2))

输出为：

ans=int((piA(1/2)*erf((1-xA2)A(1/2)))/exp(xA2),x=-1..1)

积分遇到了困难。

(选学内容一一极坐标)解决方法:

改用极坐标,也用手工确定积分限,输入:

symsrs

f=exp(-(r^2))★r;

int(int(ffrfOfl)fs,0,2*pi)

输出为:

ans=-pi★(1/exp(1)-1)

实验十无穷级数

符号表达式求和函数

symsum{S{ky)%返回符号表达式S(左)中的符号变量左从0到左-1的和值

symsum(S(k),v)%返回符号表达式S(左)中指定的符号变量左由v代替，再

对v从0至―-1的和值

symsum(S(k),a,b)%返回符号表达式S(左)中的符号变量左从a到。的和

值

symsum(S(^),v,a,b)%返回符号表达式S伏)中的符号变量左由v代替，再

对从。到。的和值

(二)符号函数的泰勒级数展开式函数

taylor(f(x))或taylor(y,x)%求函数y=/(x)的5阶麦克劳林展开式

taylor(f(x),n)%求函数/(%)的〃-1阶麦克劳林展开式

taylor(f(x),n,a)%求函数/(%)在％=a处的〃-1阶泰勒展开式

(三)在符号表达式中进行符号替换的函数

subs(S,old,new)%将符号表达式S中的符号变量old用new代替

(四)符号表达式的化简函数

simple(expr)或simplify(expr)%用于化简符号表达式expr

一、实验内容

(-)级数求和

当符号变量的和存在时，可以用symw/m命令来求无穷级数的和。

001

例1求—。

念+8〃+3

输入：

symsn

sl=symsum(1/(4*nA2+8*n+3),n,1,inf)

得到该级数的和为：

si=1/6

(二)求寨级数的收敛域

例2求£4"'白:3)"的收敛域与和函数。

〃=0〃+1

输入：

symsnx;

al=4A(2*n)*(x-3)An/(n+1);

a2=subs(al,n,n+1);

a=simplify(a2/al)

p=limit(a,n,inf)

输出为:

p=16*x-48

注意，这里对的和田都没有加绝对值。因此上式的绝对值小于1时，幕级数

收敛，大于1时发散。求出收敛区间的端点，输入：

Xl=solve(*48-x*l6=11)

X2=solve(,x*16-48=l,)

输出为：

XI=47/16

X2=49/16

由止匕可知竺<x(竺时收敛，X<”或x>竺时发散。为了判断端点的敛散

16161616

性，输入：

simplify(subs(al,'x',49/16))

得到X为右端点时塞级数的一般项为：

ans=1/(n+1)

由于级数£一~7发散，故当冗=”时发散。再输入：

+l16

simplify(subs(al,*x1,47/16))

输出为：

ans=(-1)An/(n+1)

由于级数El收敛，因此当X=%时，级数收敛。

也可以在收敛域内求得这个级数的和函数。输入：

symsnx

s3=symsum(4A(2*n)*(x-3)An/(n+l),n,0,inf)

输出为：

piecewise([49/16<=x.Inf]z[x<>49/16andabs(16*x-48)<=

1,-log(49-16*x)/(16*x-48)])

(三)将函数展开为嘉级数

例4求arctanx的5阶麦克劳林