Matlab求解微分方程组及偏微分方程组.docx
《Matlab求解微分方程组及偏微分方程组.docx》由会员分享，可在线阅读，更多相关《Matlab求解微分方程组及偏微分方程组.docx（13页珍藏版）》请在冰豆网上搜索。
Matlab求解微分方程组及偏微分方程组
第四讲Matlab求解微分方程（组）
理论介绍：
Matlab求解微分方程（组）命令
求解实例：
Matlab求解微分方程（组）实例
实际应用问题通过数学建模所归纳得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少•另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）•这就要求我们必须研究微分方程（组）的解法：
解析解法和数值解法•一•相关函数、命令及简介
1.在Matlab中，用大写字母D表示导数，Dy表示y关于自变量的一阶导数，D2y表示y关于自变量的二阶导数，依此类推.函数dsolve用来解决常微分方程
（组）的求解问题，调用格式为：
X=dsolve（‘eqn1','eqn2',…）
函数dsolve用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解.
注意，系统缺省的自变量为t
2.函数dsolve求解的是常微分方程的精确解法，也称为常微分方程的符号解.
但是，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出
其解析解，此时，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB具有丰富的函数，我们将其统称为solver，其一般格式为：
[T,Y]=solver（odefun,tspan,yO）
说明：
（1）solver为命令ode45、ode23ode113ode15sode23sode23t、ode23tb、ode15i之一.
（2）odefun是显示微分方程讨=f（t,y）在积分区间tspan=[t°,tf]上从t0到tf用初始条件y求解.
（3）如果要获得微分方程问题在其他指定时间点仏『2,川,tf上的解，则令
tspan=[t。
t1,t2J|ltf]（要求是单调的）.
（4）因为没有一种算法可以有效的解决所有的ODE问题，为此，Matlab提供
了多种求解器solver，对于不同的ODE问题，采用不同的solver.
表1Matlab中文本文件读写函数
求解器
ODE类型
特点
说明
ode45
非刚性
单步算法：
4、5阶Runge-Kutta
方程；累计截断误差（Ax）3
大部分场合的首选
算法
ode23
非刚性
单步算法：
2、3阶Runge-Kutta
方程；累计截断误差（Ax）3
使用于精度较低的
情形
ode113
非刚性
多步法：
Adams算法；高低精度
可达10“〜10》
计算时间比ode45短
ode23t
适度刚性
采用梯形算法
适度刚性情形
ode15s
刚性
多步法：
Gear's反向数值微分；
精度中等
若ode45失效时，可
尝试使用
ode23s
刚性
单步法：
2阶Rosebrock算法；低
精度
当精度较低时，计算
时间比ode15s短
ode23tb
刚性
梯形算法；低精度
当精度较低时，计算
时间比ode15s短
说明：
ode23、ode45是极其常用的用来求解非刚性的标准形式的一阶微分方
程（组）的初值问题的解的Matlab常用程序，其中：
ode23采用龙格-库塔2阶算法，用3阶公式作误差估计来调节步长，具有低等的精度•
ode45则采用龙格-库塔4阶算法，用5阶公式作误差估计来调节步长，具有中等的精度•
3.在matlab命令窗口、程序或函数中创建局部函数时，可用内联函数inline，inline函数形式相当于编写M函数文件，但不需编写M-文件就可以描述出某种数学关系调用inline函数，只能由一个matlab表达式组成，并且只能返回一个变量，不允许［u,v］这种向量形式.因而，任何要求逻辑运算或乘法运算以求得最终结果的场合，都不能应用inline函数，inline函数的一般形式为：
FunctionName=inline（函数内容'所有自变量列表'
例如：
（求解F（x）=x^2*cos（a*x）-b,a,b是标量；x是向量）在命令窗口输入：
Fofx=inline（‘x.A2*cos（a（*X）,‘x','a,'b'）;
g=Fofx（[pi/3pi/3.5],4,1）
系统输出为：
g=-1.5483-1.7259
注意：
由于使用内联对象函数inline不需要另外建立m文件，所有使用比较方便，另外在使用ode45函数的时候，定义函数往往需要编辑一个m文件来单独定义，这样不便于管理文件，这里可以使用inline来定义函数.
二•实例介绍
1.几个可以直接用Matlab求微分方程精确解的实例
例1求解微分方程y，2xy=xe'
程序：
symsxy;y=dsolve（Dy+2*x*y=x*exp（-xA2）
例2求微分方程xy'•y-ex=0在初始条件y⑴=2e下的特解并画出解函数
的图形.
程序：
symsxy;y=dsolve（x*Dy+y-exp
（1）=0'y
（1）=2*exp
（1）'''ezplot（y）
dxt
5xy=e
例3求解微分方程组
（dt在初始条件x|y=1,y|t曲=0下的特解
鱼-x—3y=0
.dt
并画出解函数的图形.
程序：
symsxyt
[x,y]=dsolve（,Dx+5*x+y=exp（t）',,Dy-x-3*y=0,,,x（0）=1,,,y（0）=0,,,t,）
simple（x）;
simple（y）
ezplot（x,y,[0,1.3]）;axisauto
2.用ode23ode45等求解非刚性标准形式的一阶微分方程（组）的初值问题
的数值解（近似解）
例4求解微分方程初值问题dx=_2y2"2x的数值解，求解范围为区
【y（0）=i
间[0,0.5].
程序：
fun=inline（,-2*y+2*xA2+2*x,,,x,,,y,）;
[x,y]=ode23（fun,[0,0.5],1）;
plot（x,y,'o-'）
例5求解微分方程
马7（1一丫2）巴y=0,y（0）=1,y'（0）=0的解，并画出
dtdt
解的图形.
分析：
这是一个二阶非线性方程，我们可以通过变换，将二阶方程化为一阶方程组求解•令Mg哼」7，则
捲（0）=1
x2（0）=0
dx!
頁％
dX2=7（1—x2）X27,.dt
编写M-文件vdp.m
functionfy=vdp（t,x）fy=[x
（2）;7*（1-x
（1）A2）*x
（2）-x
（1）];end
在Matlab命令窗口编写程序
y0=[1;0]
[t,x]=ode45（@vdp,[0,40],y0）;或[t,x]=ode45（'vdp',[0,40],y0）;y=x（:
1）;dy=x（:
2）;
plot（t,y,t,dy）
练习与思考：
M-文件vdp.m改写成inline函数程序？
3•用Euler折线法求解
Euler折线法求解的基本思想是将微分方程初值问题
f（x,y）
.y（xO=y。
化成一个代数（差分）方程，主要步骤是用差商y（xh）-y（x）替代微商鱼，于是hdx
乂二y（x。
）
记Xk^Xkh,y^y（Xk）,从而yk^y（Xkh）,于是
y。
=y（x°）,
/x“=Xk+h,k=0,1,2j||,n—1
yk1二ykhf（Xk,yQ.
例6用Euler折线法求解微分方程初值问题
dy2x
y2
dxy2
y（0）=1
的数值解（步长h取0.4）,求解范围为区间[0,2].
分析：
本问题的差分方程为
k=0,1,2,川，n-1
=0,y0=1,h=0.4
xk1=Xkh,
yk1二ykhf（Xk,yj
程序：
>>clear
>>f=sym（'y+2*x/yA2'）;
>>a=0;
>>b=2;
>>h=0.4;
>>n=（b-a）/h+1;
>>x=0;
>>y=1;
>>szj=[x,y];%数值解
>>fori=1:
n-1
y=y+h*subs（f,{'x','y'},{x,y}）;%subs，替换函数
x=x+h;
szj=[szj;x,y];
end
>>szj
>>plot（szj（:
1）,szj（:
2））
说明：
替换函数subs例如：
输入subs（a+b,a,4）意思就是把a用4替换掉,返回4+b，也可以替换多个变量，例如：
subs（cos（a）+sin（b）,{a,b},[sym（'alpha'）,2]）
分别用字符alpha替换a和2替换b,返回cos（alpha）+sin
（2）
特别说明：
本问题可进一步利用四阶Runge-Kutta法求解，Euler折线法实际
上就是一阶Runge-Kutta法，Runge-Kutta法的迭代公式为
y。
=y（x°）.
Xk厂Xkh.
h
yk1=yk■—（Li■2L2■2L3■L4）,
6
k=0,1,2,||（,n-1
Li二f（Xk,yk）,
hh
L2-f（Xk-,yk-Li）,
hh
L^f（Xk-,yk-L2）,
L4=f（Xkh,ykhL3）.
相应的Matlab程序为：
>>clear
>>f=sym（'y+2*X/yA2'）;
>>a=0;
>>b=2;
>>h=0.4;
>>n=（b-a）/h+1;
>>x=0;
>>y=1;
>>szj=[x,y];%数值解
>>fori=1:
n-1
l1=subs（f,{'x','y'},{x,y}）;替换函数
l2=subs（f,{'x','y'},{x+h/2,y+l1*h/2}）;
l3=subs（f,{'x','y'},{x+h/2,y+l2*h/2}）;
l4=subs（f,{'x','y'},{x+h,y+l3*h}）;
y=y+h*（l1+2*l2+2*l3+l4）/6;
x=x+h;
szj=[szj;x,y];
end
>>szj
>>plot（szj（:
1）,szj（:
2））
练习与思考：
（1）ode45求解问题并比较差异.
⑵利用Matlab求微分方程y⑷-2y⑶•y"=0的解.
HQI
（3）求解微分方程y—2（1-y）y•y=0,0乞x岂30,y（0）hy（0）=0的特解.
,QHII
（4）利用Matlab求微分方程初值问题（1x）y2xy,y卜卫=1,y1x^=3的解.
提醒：
尽可能多的考虑解法三•微分方程转换为一阶显式微分方程组
Matlab微分方程解算器只能求解标准形式的一阶显式微分方程（组）问题，因此在使用ODE解算器之前，我们需要做的第一步，也是最重要的一步就是借助状态变量将微分方程（组）化成Matlab可接受的标准形式.当然，如果ODEs
由一个或多个高阶微分方程给出，则我们应先将它变换成一阶显式常微分方程组下面我们以两个高阶微分方程组构成的ODEs为例介绍如何将它变换成一个一
阶显式微分方程组.
Step1将微分方程的最高阶变量移到等式左边，其它移到右边，并按阶次从低到高排列.形式为：
x（m）=f（t,x,x,x",川,x（m“y,y',y"J||,y（na））
yn）=g（t,x,x,x,川’xZyyyMyg）
升级会员 