matlab求解常微分方程docx.docx
- 文档编号:285109
- 上传时间:2022-10-08
- 格式:DOCX
- 页数:9
- 大小:49.84KB
matlab求解常微分方程docx.docx
《matlab求解常微分方程docx.docx》由会员分享,可在线阅读,更多相关《matlab求解常微分方程docx.docx(9页珍藏版)》请在冰豆网上搜索。
matlab求解常微分方程docx
用matlab求解常微分方程
在MATLAB中,由函数dsolve()解决常微分方程(组)的求解问题,其具体格式如下:
r二dsolve('eql,eq2,・••字condl,cond2,・.・;V)
匕ql,eq2,・・・*为微分方程或微分方程组,,condl,cond2,.・・;是初始条件或边界条件,P是独立变量,默认的独立变量是讥
函数dsolve用来解符号常微分方程、方程组,如果没有初始条件,则求出通解,如果有初始条件,则求出特解。
dy_1
例1:
求解常微分方程莎一石的MATLAB程序为:
dsolve(*Dy=l/(x+y)1r!
x1),
注意,系统缺省的自变量为t,因此这里要把自变量写明。
其中:
Y=lambertw(X)表示函数关系Y*exp(Y)二X。
例2:
求解常微分方程E'-y—0的MATLAB程序为:
Y2=dsolve(1y*D2y-DyA2=01,1xf)
Y2=dsolve(!
D2y*y-DyA2=0J)
我们看到有两个解,其中一个是常数0。
dx心?
—+5x+y=edt
空_兀_3『=g2f
例3:
求常微分方程组〔力'通解的MATLAB程序为:
[X,Y]=dsolve(fDx+5*x+y=exp(t),Dy-x-3*y=exp(2*t)1,111)
[X,Y]=dsolve(1Dx+2*x-Dy=l0*cos(t),Dx+Dy+2*y=4*exp(-
2*t)T,fx(0)=2,y(0)=0f,ftT)
以上这些都是常微分方程的精确解法,也称为常微分方程的符号解。
但是,我们知道,有大量的常微分方程虽然从理论上讲,其解是存在的,但我们却无法求出其解析解,此吋,我们需要寻求方程的数值解,在求常微分方程数值解方面,MATLAB具有丰富的函数,我们将其统称为solver,其一般格式为:
[T,Y]=solver(odefun,tspan,yO)
该函数表示在区间tspan=[tO,tf]±,用初始条件yO求解显式常微分方程卩=",刃。
solver为命令odc45,odc23,ode113,ode15s,ode23s,ode23t,ode23tb之一,这些命令各有特点。
我们列表说明如下:
求解器
特点
说明
ode45
一步算法,4,5阶Runge-
Kutta
方法累积截断误差(从尸
大部分场合的首选算法
ode23
一步算法,2,3阶Runge-
Kutta
方法累积截断误差(心)'
使用于精度较低的情形
odel13
多步法,Adams算法,高低精度均可达到10』~10"6
计算时间比ode45
短
ode23t
采用梯形算法
适度刚性情形
odel5s
多步法,Gear9s反向数值积分,精度中等
若ode45失效时,可尝试使用
ode23s
一步法,2阶Rosebrock算法,
低精度。
当精度较低时,计算时间比ode15s短
odefun为显式常微分方程),”,刃中的“丿)
tspan为求解区间,要获得问题在其他指定点心,也,…上的解,则令切期二[心,/“,…心](要求-单调递增或递减),yO初始条件。
例5:
求解常微分方程y'=-2y+2x2+2x,0 y=dsolve(1Dy=-2*y+2*xA2+2*x1,1y(0)=11,1x! )x=0: 0.01: 0.5; yy=subs(y,x); fun=inline(*- 2*y+2*x*x+2*x1);[x,y]=odel5s(fun,[0: 0.01: 0.5],1);ys=x・*x+exp(一 2*x); Plot(x,y,frf,x,ys,fbf) 例6: 求解常微分方程分冷+尸°"3)"'八°2()的解,并画出解的图形。 分析: 这是一个二阶非线性方程(函数以及所有偏导数军委一次幕的是现性方程,高于一次的为非线性方程),用现成的方法均不能求解,但我们可以通过下面的变换,将二阶方程化为一阶方程组,即可求解。 =dy_ 令: 兀产儿氐-亦,“=7,则得到: ~~=兀2,兀|(0)=1 at < =7(1-%,2)x2-Xj,x2(0)=0 dt 解: function[dfy]=mytt(t,fy) %fl=y;f2=dy/dt %求二阶非线性微分方程时,把一阶、二阶直到(ml)阶导数用另外一个函数代替 %用ode45命令时,必须表示成Y』f(t,Y)的形式 %Y=[yl;y2;y3],Y,=[yl\y2\y3>[y2;y3;f(yl,y2,y3)], %其中yl=y,y2=y',y3=y" %更高阶时类似 dfy=[fy (2);7*(l・fy⑴八2)*fy (2)・fy⑴]; clear;clc [t,yy]=ode45Cmytt\[040],[1;0]); plot(t,yy) legend('yTdy') 【例4.14.2.1-1]采用ODE解算扌旨令研究围绕地球旋转的卫星轨道。 (1)问题的形成 轨道上的卫星,在牛顿第二定律F=ma=m^,和万有引力定律FT叫;作用下有dtr 2 a^L=-G^r,引力常数G=6.672*10H(N.m2/kg2),ME=5.97*1024(kg)是地球的质量。 假drr 定卫星以初速度vy(0)=4000m/s在x(0)=4.2*l()7(m)处进入轨道。 (2)构成一阶微分方程组 -gme —GMe 令Y二[刃力力为]「二[尤VVxvy]T=[xy*y]T 儿 (…严 儿 (宀)计 (3)根据上式 [dYdt.m] functionYcl=DYdt(t,Y) %t %Y globalGME% xy二Y(l: 2);Vxy=Y(3: 4);% r=sqrt(sum(xy.八2)); Yd二[Vxy;-G*ME*xy//3];% (4) globalGME%<1> G=6・672e-ll;ME=5・97e24;vy0=4000;x0=-4・2e7;t0=0;tf=60*60*24*9; tspan=[tO,tf];% Y0=[xO;0;0;vyO];% X=YY(: 1); Y=YY(: 2);% plot(X,Y,,linewidth1,2);holdon %axis(*image1)% [XE,YE,ZE]=sphere(10);% RE=0・64e7;% XE=RE*XE;YE=RE*YE;ZE=0*ZE;% mesh(XE,YE,ZE),holdoff% 练习: dy_ 1.利用MATLAB求常微分方程的初值问题杰*",>? L=o=2的解。 r=dsolve(1Dy+3*y=01,! y(0)=21,1x*) 2.利用MATLAB求常微分方程的初值问题(1+扌)*=23,儿=。 二1,儿乜”的解。 r=dsolve(1D2y*(l+xA2)-2*x*Dy=0',1y(0)=1,Dy(0)=31,*) 3. 利用MATLAB求常微分方程y⑷-2厂+*=0的解。 解: y=dsolve(,D4y-2*D3y+D2y,/x,) [X,Y]=dsolve(1Dx*2+4*x+Dy- y=exp(t),Dx+3*x+y=01,1x(0)=1.5,y(0)=01,丫t】) 5.求解常微分方程/,-2(l-/)/+y=0,0 数的曲线图。 r=dsolve(1D2y-2*(l-yA2)*Dy+y=01,1y(0)=0,Dy(0)=01,*x1) functionDY=mytt2(t,Y) DY=[Y (2);2*(1-Y(l)A2)*Y (2)+Y (1)]; clear;clc [t,YY]=ode45(,mytt2\[030],[l;0]); plot(t,YY) legendC^Vdy1) 求解特殊函数方程 勒让德方程的求解 dy —2—+/(/+l)y=0dx r=dsolve(1(l-xA2)*D2y-2*x*Dy+y*l*(1+1)=01,1x') 连带勒让徳方程的求解: r=dsolve(1(l-xA2)*D2y-2*x*Dy+y*(1*(l+l)-mA2/(l-xA2))=0\) 贝塞尔方程 r=dsolve(! xA2*D2y+x*Dy+(xA2-vA2)*y=01) 利用maple mapledsolve((l-xA2)*diff(y(x),x$2)-2*x*diff(y(x)rx)+n*(n+1)*y(x)=0,y(x)) 利用MAPLE的深层符号计算资源 经典特殊函数的调用 MAPLE库函数在线帮助的检索树 发挥MAPLE的计算潜力 调用MAPLE函数 【例5.7.3.1-1]求递推方程f(n)=-3/(n-1)-2/(n-2)的通解。 (1) gsl=maple(! rsolve(f(n)=-3*f(n-1)-2*f(n-2),f(k));1) (2)调用格式二 gs2=maple(Trsolve1,Tf(n)=-3*f(n-1)-2*f(n-2)T,! f(k)! ) 【例5.7.3.1-2]求/*=供的Hessian矩阵。 (1) FHl=maple(! hessian(x*y*zf[x,y,z]);*) FH2=maple(*hessian! *x*y*z1,1[x,y,z]1) FH=sym(FH2) 【例5.7.3」・3】求sin(>2+y2)在兀二Qy=0处展开的截断8阶小量的台劳近似式。 (1) TLl=maple(^taylor(sin(xA2+yA2),[x=0,y=0],8)1) (2) maple(1readlib(mtaylor);1); TL2=maple(1mtaylor(sin(xA2+yA2),[x=0,y=0],8)1); pretty(sym(TL2)) 运行MAPLE程序 【例5.7.321】目标: 设计求取一般隐函数f(x.y)=0的导数才⑴解析解的程序,并要 求: 该程序能象MAPLE原有函数一样可被永久调用。 [DYDZZY.src] DYDZZY: =proc(f) #DYDZZY(f)isusedtogetthederivateof #animplicitfunction localEq,deq,imderiv; Eq: =1Eq1; Eq: =f; deq: =diff(Eq,x); readlib(isolate); imderiv: =isolate(deq,diff(y(x),x)); end; procread(1DYDZZY•src1) (3) sl=maple(1DYDZZY(x=log(x+y(x)));1) s2=maple(fDYDZZY(xA2*y(x)-exp(2*x)=sin(y(x)))! )s3=maple(1DYDZZY1,1cos(x+sin(y(x)))=sin(y(x))1) clea
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- matlab 求解 微分方程 docx