Matlab求解微分方程组及偏微分方程组.docx
- 文档编号:1788473
- 上传时间:2022-10-24
- 格式:DOCX
- 页数:11
- 大小:171.30KB
Matlab求解微分方程组及偏微分方程组.docx
《Matlab求解微分方程组及偏微分方程组.docx》由会员分享,可在线阅读,更多相关《Matlab求解微分方程组及偏微分方程组.docx(11页珍藏版)》请在冰豆网上搜索。
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,y0)
说明:
(1)solver为命令ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb、ode15i之一.
(2)odefun是显示微分方程在积分区间tspan上从到用初始条件求解.
(3)如果要获得微分方程问题在其他指定时间点上的解,则令tspan(要求是单调的).
(4)因为没有一种算法可以有效的解决所有的ODE问题,为此,Matlab提供了多种求解器solver,对于不同的ODE问题,采用不同的solver.
表1Matlab中文本文件读写函数
求解器
ODE类型
特点
说明
ode45
非刚性
单步算法:
4、5阶Runge-Kutta方程;累计截断误差
大部分场合的首选算法
ode23
非刚性
单步算法:
2、3阶Runge-Kutta方程;累计截断误差
使用于精度较低的情形
ode113
非刚性
多步法:
Adams算法;高低精度可达
计算时间比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.^2*cos(a*x)-b’,‘x’,’a’,’b’);
g=Fofx([pi/3pi/3.5],4,1)
系统输出为:
g=-1.5483-1.7259
注意:
由于使用内联对象函数inline不需要另外建立m文件,所有使用比较方便,另外在使用ode45函数的时候,定义函数往往需要编辑一个m文件来单独定义,这样不便于管理文件,这里可以使用inline来定义函数.
二.实例介绍
1.几个可以直接用Matlab求微分方程精确解的实例
例1求解微分方程
程序:
symsxy;y=dsolve(‘Dy+2*x*y=x*exp(-x^2)’,’x’)
例2求微分方程在初始条件下的特解并画出解函数的图形.
程序:
symsxy;y=dsolve(‘x*Dy+y-exp
(1)=0’,’y
(1)=2*exp
(1)’,’x’);ezplot(y)
例3求解微分方程组在初始条件下的特解并画出解函数的图形.
程序:
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.用ode23、ode45等求解非刚性标准形式的一阶微分方程(组)的初值问题的数值解(近似解)
例4求解微分方程初值问题的数值解,求解范围为区间[0,0.5].
程序:
fun=inline('-2*y+2*x^2+2*x','x','y');
[x,y]=ode23(fun,[0,0.5],1);
plot(x,y,'o-')
例5求解微分方程的解,并画出解的图形.
分析:
这是一个二阶非线性方程,我们可以通过变换,将二阶方程化为一阶方程组求解.令,则
编写M-文件vdp.m
functionfy=vdp(t,x)
fy=[x
(2);7*(1-x
(1)^2)*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折线法求解的基本思想是将微分方程初值问题
化成一个代数(差分)方程,主要步骤是用差商替代微商,于是
记从而于是
例6用Euler折线法求解微分方程初值问题
的数值解(步长取0.4),求解范围为区间[0,2].
分析:
本问题的差分方程为
程序:
>>clear
>>f=sym('y+2*x/y^2');
>>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法的迭代公式为
相应的Matlab程序为:
>>clear
>>f=sym('y+2*x/y^2');
>>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求解问题并比较差异.
(2)利用Matlab求微分方程的解.
(3)求解微分方程的特解.
(4)利用Matlab求微分方程初值问题的解.
提醒:
尽可能多的考虑解法
三.微分方程转换为一阶显式微分方程组
Matlab微分方程解算器只能求解标准形式的一阶显式微分方程(组)问题,因此在使用ODE解算器之前,我们需要做的第一步,也是最重要的一步就是借助状态变量将微分方程(组)化成Matlab可接受的标准形式.当然,如果ODEs由一个或多个高阶微分方程给出,则我们应先将它变换成一阶显式常微分方程组.下面我们以两个高阶微分方程组构成的ODEs为例介绍如何将它变换成一个一阶显式微分方程组.
Step1将微分方程的最高阶变量移到等式左边,其它移到右边,并按阶次从低到高排列.形式为:
Step2为每一阶微分式选择状态变量,最高阶除外
注意:
ODEs中所有是因变量的最高阶次之和就是需要的状态变量的个数,最高阶的微分式不需要给它状态变量.
Step3根据选用的状态变量,写出所有状态变量的一阶微分表达式
练习与思考:
(1)求解微分方程组
其中
(2)求解隐式微分方程组
提示:
使用符号计算函数solve求,然后利用求解微分方程的方法
四.偏微分方程解法
Matlab提供了两种方法解决PDE问题,一是使用pdepe函数,它可以求解一般的PDEs,具有较大的通用性,但只支持命令形式调用;二是使用PDE工具箱,可以求解特殊PDE问题,PDEtoll有较大的局限性,比如只能求解二阶PDE问题,并且不能解决片微分方程组,但是它提供了GUI界面,从复杂的编程中解脱出来,同时还可以通过File—>SaveAs直接生成M代码.
1.一般偏微分方程(组)的求解
(1)Matlab提供的pdepe函数,可以直接求解一般偏微分方程(组),它的调用格式为:
sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)
@pdefun是PDE的问题描述函数,它必须换成标准形式:
这样,PDE就可以编写入口函数:
[c,f,s]=pdefun(x,t,u,du),m,x,t对应于式中相关参数,du是u的一阶导数,由给定的输入变量可表示出c,f,s这三个函数.
@pdebc是PDE的边界条件描述函数,它必须化为形式:
于是边值条件可以编写函数描述为:
[pa,qa,pb,qb]=pdebc(x,t,u,du),其中a表示下边界,b表示上边界.
@pdeic是PDE的初值条件,必须化为形式:
,故可以使用函数描述为:
u0=pdeic(x)
sol是一个三维数组,sol(:
:
i)表示的解,换句话说,对应x(i)和t(j)时的解为sol(i,j,k),通过sol,我们可以使用pdeval函数直接计算某个点的函数值.
(2)实例说明
求解偏微分
其中,且满足初始条件及边界条件
解:
(1)对照给出的偏微分方程和pdepe函数求解
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Matlab 求解 微分 方程组