matlab求解常微分方程docxWord下载.docx
《matlab求解常微分方程docxWord下载.docx》由会员分享，可在线阅读，更多相关《matlab求解常微分方程docxWord下载.docx（9页珍藏版）》请在冰豆网上搜索。
例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-
方法累积截断误差（心）'
使用于精度较低的情形
odel13
多步法，Adams算法，高低精度均可达到10』~10"
6
计算时间比ode45
短
ode23t
采用梯形算法
适度刚性情形
odel5s
多步法，Gear9s反向数值积分，精度中等
若ode45失效时,可尝试使用
ode23s
一步法，2阶Rosebrock算法，
低精度。
当精度较低时，计算时间比ode15s短
odefun为显式常微分方程）,”，刃中的“丿）
tspan为求解区间，要获得问题在其他指定点心，也，…上的解，则令切期二［心,/“，…心］（要求-单调递增或递减），yO初始条件。
例5：
求解常微分方程y'
=-2y+2x2+2x,0<
x<
0.5,y（0）=l的MATLAB程序如下：
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.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;
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<
30,y（0）=l,/（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）];
[t,YY]=ode45（,mytt2\[030],[l;
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矩阵。
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阶小量的台劳近似式。
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;
=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
升级会员 