Matlab求解微分方程组及偏微分方程组.docx
- 文档编号:29192181
- 上传时间:2023-07-21
- 格式:DOCX
- 页数:15
- 大小:73.41KB
Matlab求解微分方程组及偏微分方程组.docx
《Matlab求解微分方程组及偏微分方程组.docx》由会员分享,可在线阅读,更多相关《Matlab求解微分方程组及偏微分方程组.docx(15页珍藏版)》请在冰豆网上搜索。
Matlab求解微分方程组及偏微分方程组
第四讲Matlab求解微分方程(组)
理论介绍:
Matlab求解微分方程(组)命令
求解实例:
Matlab求解微分方程(组)实例
实际应用问题通过数学建模所归纳得到的方程,绝大多数都是微分方程,真正能得到代数方程的机会很少.另一方面,能够求解的微分方程也是十分有限的,特别是高阶方程和偏微分方程(组).这就要求我们必须研究微分方程(组)的解法:
解析解法和数值解法.
一.相关函数、命令及简介
1.在Matlab中,用大写字母D表示导数,Dy表示y关于自变量的一阶导数,D2y表示y关于自变量的二阶导数,依此类推.函数dsolve用来解决常微分方程
(组)的求解问题,调用格式为:
X=dsolve(<eqnl,,,eqn2
函数dsolve用来解符号常微分方程、方程组,如果没有初始条件,则求出通解,如果有初始条件,则求出特解.
注意,系统缺省的自变量为t
2.函数dsolve求解的是常微分方程的精确解法,也称为常微分方程的符号解.但是,有大量的常微分方程虽然从理论上讲,其解是存在的,但我们却无法求出其解析解,此时,我们需要寻求方程的数值解,在求常微分方程数值解方面,MATLAB具有丰富的函数,我们将其统称为solver,其一般格式为:
[T,Y]=solver(odefun,tspan,yO)
说明:
(1)solver为命令ode45、ode23、odel13、odel5s、ode23s、ode23t、ode23tb、odel5i之一.
(2)odefun是显示微分方程),=f(t,y)在积分区间tspan=[心心]上从心到“用初始条件儿求解.
(3)如果要获得微分方程问题在其他指定时间点bG©…心上的解,则令
(span=『“,•••『/■](要单调的).
(4)因为没有一种算法可以有效的解决所有的ODE问题,为此,Matlab提供T多种求解器solver,对于不同的ODE问题,采用不同的solver.
表1Matlab中文本文件读写函数
求解器
ODE类型
特点
说明
ode45
非刚性
单步算法:
4、5阶Runge-Kutta
方程;累计截断误差(Ax)'
大部分场合的首选
算法
ode23
非刚性
单步算法:
2、3阶Runge-Kutta
方程;累计截断误差(心)'
使用于精度较低的
情形
odell3
非刚性
多步法:
Adams算法;高低精度
可达ICT?
~10"
计算时间比ode45短
ode23t
适度刚性
采用梯形算法
适度刚性情形
ode15s
刚性
多步法:
Gear's反向数值微分;
精度中等
若ode45失效时,可
尝试使用
ode23s
刚性
单步法:
2阶Rosebrock算法;低
精度
当精度较低时,计算
时间比ode15s短
ode23tb
刚性
梯形算法;低精度
当精度较低时,计算
时间比ode15s短
说明:
。
血23、。
血45是极其常用的用来求解非刚性的标准形式的一阶微分方
程(组)的初值问题的解的Matlab常用程序,其中:
ode23采用龙格-库塔2阶算法,用3阶公式作误差估计来调节步长,具有低等的精度.
。
血45则采用龙格-库塔4阶算法,用5阶公式作误差估计来调节步长,具有中等的精度.
3.在matlab命令窗口、程序或函数中创建局部函数时,可用联函数inline,inline函数形式相当于编写M函数文件,但不需编写文件就可以描述出某种数学关系.调用inline函数,只能由一个matlab表达式组成,并且只能返回一个变量,不允许[u,v]这种向量形式.因而,任何要求逻辑运算或乘法运算以求得最终结果的场合,都不能应用inline函数,inline函数的一般形式为:
FunctionName=inline(*函数容','所有自变量列表')
例如:
(求解F(x)=xA2*cos(a*x)-b,a,b是标量;x是向量)在命令窗口输入:
Fofk=inline(6x・A2*cos(a*x)-b>,仪‘七丁!
/);
g=Fofx([pi/3pi/3・5],4J)
系统输出为:
g=-1.5483-1.7259
注意:
山于使用联对象函数inline不需要另外建立m文件,所有使用比较方便,另外在使用ode45函数的时候,定义函数往往需要编辑一个m文件来单独定义,这样不便于管理文件,这里可以使用inline来定义函数.
实例介绍
1.儿个可以直接用Matlab求微分方程精确解的实例
例1求解微分方程y+2卩,=
程序:
symsxy;y=dsolve(tDy+2*x*y=x*exp(-xA2)\,x,)
例2求微分方程厂0在初始条件y(\)=2e下的特解并画出解函数的图形.
程序:
symsxy;y=dsolve('x*Dy+y・exp(l)=0','y(l)=2*exp(l)','x');ezpk>t(y)
dx,
—+5x+y=e
例3求解微分方程组山~在初始条件H『l,ylj=O下的特解^-x-3y=0
[dt
并画出解函数的图形.
程序:
symsxyt
[x,y]=dsolve(,Dx+5*x+y=exp(t),/Dy-x-3*y=0\'x(0)=lVy(0)=0Vt)
simple(x);
simple(y)
ezplot(x,y,[0,1.3]);axisauto
2.用ode23.ode45等求解非刚性标准形式的一阶微分方程(组)的初值问题的数值解(近似解)
(小,
丿=v-I-?
V24-?
V
例4求解微分方程初值问题页一>的数值解,求解围为区间
y(0)=1
[0,0.5].
程序:
fun=inline('-2*y+2*xA2+2*x','x',,y,);
[x,y]=ode23(fun,[0,0.5],1);plot(x,y/oJ)
例5求解微分方程[律一"(I")牛+y=0y(0)=],y(0)=0的解,并画出drdt
解的图形.
分析:
这是一个二阶非线性方程,我们可以通过变换,将二阶方程化为一阶
方程组求解•令x1=y,x2=
—,//=7,则dt
字=禺,X,(0)=1
dt・
dx
—=7(1-X;)x7一X],as(0)=0dt
编写M■文件vdp.m
functionfy=vdp(t,x)
fy=[x
(2);7*(1-x(l)A2)*x⑵-x
(1)];
end
在Matlab命令窗「I编写程序
yo=[1;oj
[t,x]=ode45(vdp,[0,40],y0);或[t,x]=ode45('vdp;[0,40],y0);
y=x(:
l);dy=x(:
2);
plot(t,y,t,dy)
练习与思考:
M-文件vdp.m改写成inline函数程序?
3.用Euler折线法求解
Euler折线法求解的基本思想是将微分方程初值问题
dy—、
.y(^o)=>o 化成一个代数(差分)方程,主要步骤是用差商+")—)'(')替代微商空,于是 hdx =yM 记札I=Xk+h,yk=y(xk),从而yk+}=y(xk+/? ),于是 儿=〉'(观), .儿+严儿+妙幺,儿)・ 例6用Euler折线法求解微分方程初值问题 dy2x —=y+— y(0)=1 的数值解(步长力取0.4),求解围为区间[0,2]. 分析: 本问题的差分方程为 Xq=0,y°=1、h=0.4 .儿+严儿+妙’(兀,儿)• 程序: »clear »f=sym(,y+2*x/yA2,); »a=0; »b=2; »h=0.4; »n=(b-a)/h+l; »x=0; »y=i; »szj=[x,y];%数值解 »fori=l: n-l y二y+h*subs(f,{'x','y'},{x,y});%subs,替换函数 x=x+h; szj二[szj;x,y]; end »szj »plot(szj(: J),szj(: 2)) 说明: 替换函数subs例如: 输入siibs(“+b,a,4)意思就是把d用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法的迭代公式为 儿=)仇), 和=耳+九 儿+i=以+匸(厶+2厶j+2厶+厶4), O 厶 厶=f(xk+2”+? 厶), k=0丄2,…,”一1 厶=/(xt+—074■—厶)' L4=f(xk+h,yk+hLi). 相应的Matlab程序为: »clear »f=sym(,y+2*x/yA2,); »a=0; »b=2; »h=0.4; »n=(b-a)/h+l; »x=0; »y=l; »szj=[x,y];%数值解 »for上l: n-l ll=subs(f,「x;y},{x,y});替换函数 12=subs(f,{'x'/y'),{x+h/2,y+l1*h/2}); 13=subs(f,{'x'/y'),{x+h/2,y+12*h/2}); 14=subs(f,{'xTy'}9{x+h、y+13*h}); y=y+h*(ll+2*12+2*13+14)/6; x=x+h; szj二[szj;x,y]; end »szj »plot(szj(: l),szj(: 2)) 练习与思考: (1)。 血45求解问题并比较差异. (2)利用Matlab求微分方程y⑷-2/3)+/=0的解. (3)求解微分方程一2(1-y2)y+y=0,0 (4)利用Matlab求微分方程初值问题(1+x2)/=2xy\yl“=1,yI“=3的解. 提醒: 尽可能多的考虑解法 三.微分方程转换为一阶显式微分方程组 Matlab微分方程解算器只能求解标准形式的一阶显式微分方程(组)问题,因此在使用ODE解算器之前,我们需要做的笫一步,也是最重要的一步就是借助状态变量将微分方程(组)化成Matlab可接受的标准形式.当然,如果ODEs山一个或多个高阶微分方程给出,则我们应先将它变换成一阶显式常微分方程组.下面我们以两个高阶微分方程组构成的ODEs为例介绍如何将它变换成一个一阶显式微分方程组. Stepl将微分方程的最高阶变量移到等式左边,其它移到右边,并按阶次从低到高排列•形式为: Step2为每一阶微分式选择状态变量,最高阶除外 注意: ODEs中所有是因变量的最高阶次之和就是需要的状态变量的个数,最高阶的微分式不需要给它状态变量. Step3根据选用的状态变量,写出所有状态变量的一阶微分表达式 片=X2,X2=召,兀3兀“=/(人召 •■ 兀1+1—•'加十2,•••,兀T—g(人X],兀29尤3,•••9) 练习与思考: (1)求解微分方程组 其中r2=Ja_”)2+y2,斤=J(x+〃)2+y2,K=1一“,//=1/82.45,x(0)=1.2,y(0)=0,x(0)=0,y(0)=-1.049355751 (2)求解隐式微分方程组 x+2yx=2y < ■■9•・ xy+3xy+xy_y=5 提示: 使用符号计算函数solve求V,y\然后利用求解微分方程的方法 四.偏微分方程解法 Matlab提供了两种方法解决PDE问题,一是使用pdepe函数,它可以求解一般的PDEs,具有较大的通用性,但只支持命令形式调用;二是使用PDE工具箱,可以求解特殊PDE问题,PDEtoll有较大的局限性,比如只能求解二阶PDE问题,并且不能解决片微分方程组,但是它提供了GUI界面,从复杂的编程中解脱出来,同时还可以通过File一>SaveAs直接生成M代码. 1.一般偏微分方程(组)的求解 (l)Matlab提供的pdepe函数,可以直接求解一般偏微分方程(组),它的调用格式为: sol=pdepe(m,pdefun,pdeic,pdebc,x,t) pdefun是PDE的问题描述函数,它必须换成标准形式: /,duxdud.,dir..,dir c(xj,—)—=x—|.vf(x,t,w,—)]+s(x,t,u,—) oxotoxoxox 这样,PDE就可以编写入口函数: [c,f,s]=pdefun(x,t,u,du),m,x,t对应于式中相关参数,du是u的一阶导数,山给定的输入变量可表示出c,f,s这三个函数. pdebc是PDE的边界条件描述函数,它必须化为形式: 6以 p(xj.u)=q(x,tyu)*f(x,t,u,—)=0 ox 于是边值条件可以编写函数描述为: [pa,qa,pb,qb]二pdebc(x,t,u,du),其中a表示下边界,b表示上边界. pdeic是PDE的初值条件,必须化为形式: h(xJ0)=//0,故可以使用函数描 述为: uO二pdeic(x) sol是一个三维数组,sol(: ;i)表示比的解,换句话说,你对应x(i)和t(j)时的 解为solgk),通过sol,我们可以使用pdeval函数直接计算某个点的函数值. (2)实例说明 求解偏微分 %□标PDE函数 function[c,f,s]=pdefun(x,t,uJu) C=[l;l]; f=[0.024*du(l);0.17*du (2)]; temp=u(l)-u (2); s=[-l;l].*(exp(5.73*temp)-exp(-l1.46*temp)) end (2)边界条件改写为: %边界条件函数 function[pa,qa,pb5qb]=pdebc(xa,ua,xb,ub,t) pa=[0;ua⑵]; qa=[l;0]; pb=[ub(l)-l;O]; qb=[O;l]; end (3)初值条件改写为: %初值条件函数 functionuO=pdeic(x) u0=[l;0]; end (4)编写主调函数 clcx=0: 0.05: l; t=0: 0.05: 2; m=0; sol=pdepe(m,pdefun,pdeic,pdebc,x,t); subplot(2JJ)surf(x,t,sol(: : J)) subplot(2,l,2)surf(x,t,sol(: : 2)) 练习与思考: Thisexampleillustratesthestraightforwardformulation,computation,andplottingofthesolutionofasinglePDE・ ,duddii 矿—=一 (一) dtdxdx Thisequationholdsonaninterval0 "(0./)=0;/rk+——(lj)=0 dx 2.PDEtool求解偏微分方程 (1)PDEtool(GUI)求解偏微分方程的一般步骤 在Matlab命令窗口输入pdetool,回车,PDE工具箱的图形用户界面(GUI)系统就启动了•从定义一个偏微分方程问题到完成解偏微分方程的定解,整个过程大致可以分为六个阶段 Step1“Draw模式”绘制平面有界区域C,通过公式把Matlab系统提供的实体模型: 矩形、圆、椭圆和多边形,组合起来,生成需要的平面区域. Step2"Boundary模式”定义边界,声明不同边界段的边界条件. Step3“PDE模式”定义偏微分方程,确定方程类型和方程系数c,a,f,d,根据具体情况,还可以在不同子区域声明不同系数. Step4"Mesh模式”网格化区域G,可以控制自动生成网格的参数,对生成的网格进行多次细化,使网格分割更细更合理. Step5“Solve模式”解偏微分方程,对于椭圆型方程可以激活并控制非线性自适应解题器来处理非线性方程;对于抛物线型方程和双曲型方程,设置初始边界条件后可以求出给定时刻t的解;对于特征值问题,可以求出给定区间上的特征值.求解完成后,可以返回到Step4,对网格进一步细化,进行再次求解. Step6“View模式”计算结果的可视化,可以通过设置系统提供的对话框,显示所求的解的表面图、网格图、等高线图和箭头梯形图.对于抛物线型和双曲线型问题的解还可以进行动画演示. (2)实例说明用法 求解一个正方形区域上的特征值问题: 一△“u=加 <2 正方形区域为: —1K1,—1K1. (1)使用PDE工具箱打开GUI求解方程 (2)进入Draw模式,绘制一个矩形,然后双击矩形,在弹出的对话框中设置 Left=-1,Bottom=-l,Width=2,Height=2,确认并关闭对话框 (3)进入Boundary模式,边界条件采用Dirichlet条件的默认值 (4)进入PDE模式,单击工具栏PDE按钮,在弹出的对话框中方程类型选择 Eigenmodes,参数设置c=1,a=-l/2,d=1,确认后关闭对话框 (5)单击工具栏的△按钮,对正方形区域进行初始网格剖分,然后再对网格进一步细化剖分一次 (6)点开solve菜单,单击Parameters选项,在弹出的对话框中设置特征值区域为卜20,20] (7)单击Plot菜单的Parameters项,在弹出的对话框中选中Color、Height(3-Dplot)和showmesh项,然后单击Done确认 (8)单击工具栏的“二”按钮,开始求解
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Matlab 求解 微分 方程组