Matlab求解微分方程(组)及偏微分方程(组).pdf

该【Matlab求解微分方程(组)及偏微分方程(组) 】是由【我是开始】上传分享，文档一共【13】页，该文档可以免费在线阅读，需要了解更多关于【Matlab求解微分方程(组)及偏微分方程(组) 】的内容，可以使用淘豆网的站内搜索功能，选择自己适合的文档，以下文字是截取该文章内的部分文字，如需要获得完整电子版，请下载此文档到您的设备，方便您编辑和打印。:.
创作编号:BG75314000**********SX
创作者:别如克*
第四讲Matlab求解微分方程(组)
理论介绍:Matlab求解微分方程(组)命令
求解实例:Matlab求解微分方程(组)实例
实际应用问题通过数学建模所归纳得到的方程,绝大多数都是微分方
程,,能够求解的微分方程也是
十分有限的,特别是高阶方程和偏微分方程(组).这就要求我们必须研究
微分方程(组)的解法:解析解法和数值解法.
、命令及简介
,用大写字母D表示导数,Dy表示y关于自变量的一阶
导数,D2y表示y关于自变量的二阶导数,
常微分方程(组)的求解问题,调用格式为:
X=dsolve(‘eqn1’,’eqn2’,…)
函数dsolve用来解符号常微分方程、方程组,如果没有初始条件,则
求出通解,如果有初始条件,则求出特解.
注意,系统缺省的自变量为t
,也称为常微分方程的符
,有大量的常微分方程虽然从理论上讲,其解是存在的,但我们
却无法求出其解析解,此时,我们需要寻求方程的数值解,在求常微分方
程数值解方面,MATLAB具有丰富的函数,我们将其统称为solver,其一
般格式为:
[T,Y]=solver(odefun,tspan,y0)
说明:(1)solver为命令ode45、ode23、ode113、ode15s、ode23s、ode23t、
ode23tb、ode15i之一.
1:.
(2)odefun是显示微分方程y'f(t,y)在积分区间tspan[t,t]上从t
0f0
到t用初始条件y求解.
f0
(3)如果要获得微分方程问题在其他指定时间点t,t,t,,t上的解,则
012f
令tspan[t,t,t,t](要求是单调的).
012f
(4)因为没有一种算法可以有效的解决所有的ODE问题,为此,Matlab
提供了多种求解器solver,对于不同的ODE问题,采用不同的solver.
表1Matlab中文本文件读写函数
求解器ODE类型特点说明
单步算法:4、5阶Runge-Kutta
大部分场合的首选
ode45非刚性
方程;累计截断误差(x)3算法
单步算法:2、3阶Runge-Kutta
使用于精度较低的
ode23非刚性
方程;累计截断误差(x)3情形
多步法:Adams算法;高低精度
ode113非刚性计算时间比ode45短
可达103~106
ode23t适度刚性采用梯形算法适度刚性情形
多步法:Gear’s反向数值微分;若ode45失效时,可
ode15s刚性
精度中等尝试使用
单步法:2阶Rosebrock算法;低当精度较低时,计算
ode23s刚性
精度时间比ode15s短
当精度较低时,计算
ode23tb刚性梯形算法;低精度
时间比ode15s短
说明:ode23、ode45是极其常用的用来求解非刚性的标准形式的一阶
微分方程(组)的初值问题的解的Matlab常用程序,其中:
ode23采用龙格-库塔2阶算法,用3阶公式作误差估计来调节步长,
具有低等的精度.
ode45则采用龙格-库塔4阶算法,用5阶公式作误差估计来调节步长,
具有中等的精度.
2:.
、程序或函数中创建局部函数时,可用内联函数
inline,inline函数形式相当于编写M函数文件,但不需编写M-文件就可以
,只能由一个matlab表达式组成,并
且只能返回一个变量,不允许[u,v],任何要求逻辑运算
或乘法运算以求得最终结果的场合,都不能应用inline函数,inline函数的
一般形式为:
FunctionName=inline(‘函数内容’,‘所有自变量列表’)
例如:(求解F(x)=x^2*cos(a*x)-b,a,b是标量;x是向量)在命令窗口输入:
Fofx=inline(‘x.^2*cos(a*x)-b’,‘x’,’a’,’b’);
g=Fofx([pi/3pi/],4,1)
系统输出为:g=--
创作编号:BG75314000**********SX
创作者:别如克*
注意:由于使用内联对象函数inline不需要另外建立m文件,所有使
用比较方便,另外在使用ode45函数的时候,定义函数往往需要编辑一个
m文件来单独定义,这样不便于管理文件,这里可以使用inline来定义函
数.


例1求解微分方程y'2xyxex2
程序:symsxy;y=dsolve(‘Dy+2*x*y=x*exp(-x^2)’,’x’)
例2求微分方程xy'yex0在初始条件y(1)2e下的特解并画出解
函数的图形.
程序:symsxy;
y=dsolve(‘x*Dy+y-exp(1)=0’,’y(1)=2*exp(1)’,’x’);ezplot(y)
3:.
dx
5xyet

dt
例3求解微分方程组在初始条件x|1,y|0下的
dyt0t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,]);axisauto
、ode45等求解非刚性标准形式的一阶微分方程(组)的初
值问题的数值解(近似解)
dy
2y2x22x
例4求解微分方程初值问题dx的数值解,求解范围
y(0)1

为区间[0,].
程序:fun=inline('-2*y+2*x^2+2*x','x','y');
[x,y]=ode23(fun,[0,],1);
plot(x,y,'o-')
d2ydy
例5求解微分方程(1y2)y0,y(0)1,y'(0)0的解,并
dt2dt
画出解的图形.
分析:这是一个二阶非线性方程,我们可以通过变换,将二阶方程化
dy
y,x,7,则
12dt
dx
1x,x(0)1
21
dt

dx
27(1x2)xx,x(0)0
dt1212
编写M-
functionfy=vdp(t,x)
fy=[x(2);7*(1-x(1)^2)*x(2)-x(1)];
4:.
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-?

Euler折线法求解的基本思想是将微分方程初值问题
dy
f(x,y)
dx
y(x)y

00
y(xh)y(x)dy
化成一个代数(差分)方程,主要步骤是用差商替代微商,
hdx
于是
y(xh)y(x)
kkf(x,y(x))
hkk
yy(x)

00
记xxh,yy(x),从而yy(xh),于是
k1kkkk1k
yy(x),
00

xxh,k0,1,2,,n1
k1k

yyhf(x,y).
k1kkk
例6用Euler折线法求解微分方程初值问题
dy2x
y
dxy2

y(0)1
的数值解(),求解范围为区间[0,2].
分析:本问题的差分方程为
x0,y1,h
00
xxh,k0,1,2,,n1
k1k

yyhf(x,y).
k1kkk
5:.
程序:>>clear
>>f=sym('y+2*x/y^2');
>>a=0;
>>b=2;
>>h=;
>>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),
00

xxh,
k1k
h
yy(L2L2LL),
k1k61234

Lf(x,y),
1kkk0,1,2,,n1
hh
Lf(x,yL),
2k2k21

hh
Lf(x,yL),
3k2k22

Lf(xh,yhL).
4kk3
6:.
相应的Matlab程序为:>>clear
>>f=sym('y+2*x/y^2');
>>a=0;
>>b=2;
创作编号:BG75314000**********SX
创作者:别如克*
>>h=;
>>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求解问题并比较差异.
(2)利用Matlab求微分方程y(4)2y(3)y''0的解.
(3)求解微分方程y''2(1y2)y'y0,0x30,y(0)1,y,(0)0的特
7:.
解.
(4)利用Matlab求微分方程初值问题(1x2)y''2xy',y|1,y'|3的
x0x0
解.
提醒:尽可能多的考虑解法

Matlab微分方程解算器只能求解标准形式的一阶显式微分方程(组)
问题,因此在使用ODE解算器之前,我们需要做的第一步,也是最重要的
一步就是借助状态变量将微分方程(组)化成Matlab可接受的标准形式.
当然,如果ODEs由一个或多个高阶微分方程给出,则我们应先将它变换

例介绍如何将它变换成一个一阶显式微分方程组.
Step1将微分方程的最高阶变量移到等式左边,其它移到右边,并按
:
x(m)f(t,x,x',x'',,x(m1),y,y',y'',,y(n1))

y(n)g(t,x,x',x'',,x(m1),y,y',y'',,y(n1))

Step2为每一阶微分式选择状态变量,最高阶除外
xx,xx',xx'',,xx(m1),
123m
xy,xy',xy'',,xy(n1)
m1m2m3mn
注意:ODEs中所有是因变量的最高阶次之和就是需要的状态变量的个
数,最高阶的微分式不需要给它状态变量.
Step3根据选用的状态变量,写出所有状态变量的一阶微分表达式
x'x,x'x,x'x,,x'f(t,x,x,x,,x)
122334m123mn
x'x,,x'g(t,x,x,x,,x)
m1m2mn123mn
练****与思考:(1)求解微分方程组
*(x)(x*)
x''2y'x

r3r3
12
*yy
'''
y2xy
r3r3

12
其中
r(x*)2y2,r(x)2y2,*1,1/,x(0),
21
8:.
y(0)0,x'(0)0,y'(0)
(2)求解隐式微分方程组
x''2y'x2y''

x''y'3x'y''xy'y5

提示:使用符号计算函数solve求x'',y'',然后利用求解微分方程的方
法

Matlab提供了两种方法解决PDE问题,一是使用pdepe函数,它可以
求解一般的PDEs,具有较大的通用性,但只支持命令形式调用;二是使用
PDE工具箱,可以求解特殊PDE问题,PDEtoll有较大的局限性,比如只
能求解二阶PDE问题,并且不能解决片微分方程组,但是它提供了GUI
界面,从复杂的编程中解脱出来,同时还可以通过File—>SaveAs直接生
成M代码.
(组)的求解
(1)Matlab提供的pdepe函数,可以直接求解一般偏微分方程(组),它
的调用格式为:sol=pdepe(m,***@pdefun,***@pdeic,***@pdebc,x,t)
***@pdefun是PDE的问题描述函数,它必须换成标准形式:
uuuu
c(x,t,)xm[xmf(x,t,u,)]s(x,t,u,)
xtxxx
这样,PDE就可以编写入口函数:[c,f,s]=pdefun(x,t,u,du),m,x,t对应于式
中相关参数,du是u的一阶导数,由给定的输入变量可表示出c,f,s这三个
函数.
***@pdebc是PDE的边界条件描述函数,它必须化为形式:
u
p(x,t,u)q(x,t,u).*f(x,t,u,)0
x
于是边值条件可以编写函数描述为:[pa,qa,pb,qb]=pdebc(x,t,u,du),其中a表
示下边界,b表示上边界.
***@pdeic是PDE的初值条件,必须化为形式:u(x,t)u,故可以使用
00
函数描述为:u0=pdeic(x)
sol是一个三维数组,sol(:,:,i)表示u的解,换句话说,u对应x(i)和t(j)
ik
9:.
时的解为sol(i,j,k),通过sol,我们可以使用pdeval函数直接计算某个点的
函数值.
(2)实例说明
求解偏微分
u2u
1F(uu)
tx212

u2u

2F(uu)
tx212
其中,F(x)e(x,0)1,u(x,0)0及边界条件
12
uu
1(0,t)0,u(0,t)0,u(1,t)1,2(1,t)0
x21x
解:(1)对照给出的偏微分方程和pdepe函数求解的标准形式,原方程
改写为
u

1uF(uu)
x
.*112

1tuxuF(uu)

x
u

1F(uu)
x
可见m0,c,f,s12

1uF(uu)

x
%目标PDE函数
function[c,f,s]=pdefun(x,t,u,du)
c=[1;1];
f=[*du(1);*du(2)];
temp=u(1)-u(2);
s=[-1;1].*(exp(*temp)-exp(-*temp))
end
(2)边界条件改写为:
010u110
下边界.*f上边界1.*f

u00000
2
%边界条件函数
10:.
function[pa,qa,pb,qb]=pdebc(xa,ua,xb,ub,t)
pa=[0;ua(2)];
qa=[1;0];
pb=[ub(1)-1;0];
qb=[0;1];
end
u1
(3)初值条件改写为:1

u0
2
%初值条件函数
functionu0=pdeic(x)
创作编号:BG75314000**********SX
创作者:别如克*
u0=[1;0];
end
(4)编写主调函数
clc
x=0::1;
t=0::2;
m=0;
sol=pdepe(m,***@pdefun,***@pdeic,***@pdebc,x,t);
subplot(2,1,1)surf(x,t,sol(:,:,1))
subplot(2,1,2)surf(x,t,sol(:,:,2))
练****与思考:Thisexampleillustratesthestraightforwardformulation,
computation,andplottingofthesolutionofasinglePDE.
uu
2()
txx
11:.
Thisequationholdsonaninterval0x1fortimest
theinitialconditionu(x,0)sinxandboundaryconditions
u
u(0,t)0;et(1,t)0
x

(1)PDEtool(GUI)求解偏微分方程的一般步骤
在Matlab命令窗口输入pdetool,回车,PDE工具箱的图形用户界面
(GUI)
解,整个过程大致可以分为六个阶段
Step1“Draw模式”绘制平面有界区域,通过公式把Matlab系统
提供的实体模型:矩形、圆、椭圆和多边形,组合起来,生成需要的平面
区域.
Step2“Boundary模式”定义边界,声明不同边界段的边界条件.
Step3“PDE模式”定义偏微分方程,确定方程类型和方程系数c,a,f,d,
根据具体情况,还可以在不同子区域声明不同系数.
Step4“Mesh模式”网格化区域,可以控制自动生成网格的参数,
对生成的网格进行多次细化,使网格分割更细更合理.
Step5“Solve模式”解偏微分方程,对于椭圆型方程可以激活并控制
非线性自适应解题器来处理非线性方程;对于抛物线型方程和双曲型方程,
设置初始边界条件后可以求出给定时刻t的解;对于特征值问题,可以求
,可以返回到Step4,对网格进一步细
化,进行再次求解.
Step6“View模式”计算结果的可视化,可以通过设置系统提供的对
话框,显示所求的解的表面图、网格图、
线型和双曲线型问题的解还可以进行动画演示.
(2)实例说明用法
求解一个正方形区域上的特征值问题:
1
uuu
2
u|0


12:.
正方形区域为:1x1,1x1.
(1)使用PDE工具箱打开GUI求解方程
(2)进入Draw模式,绘制一个矩形,然后双击矩形,在弹出的对话框
中设置Left=-1,Bottom=-1,Width=2,Height=2,确认并关闭对话框
(3)进入Boundary模式,边界条件采用Dirichlet条件的默认值
(4)进入PDE模式,单击工具栏PDE按钮,在弹出的对话框中方程类
型选择Eigenmodes,参数设置c=1,a=-1/2,d=1,确认后关闭对话框
(5)单击工具栏的按钮,对正方形区域进行初始网格剖分,然后再对
网格进一步细化剖分一次
(6)点开solve菜单,单击Parameters选项,在弹出的对话框中设置特征
值区域为[-20,20]
(7)单击Plot菜单的Parameters项,在弹出的对话框中选中Color、
Height(3-Dplot)和showmesh项,然后单击Done确认
(8)单击工具栏的“=”按钮,开始求解
创作编号:BG75314000**********SX
创作者:别如克*
13