Matlab求解微分方程组及偏微分方程组.docx

1、Matlab求解微分方程组及偏微分方程组第四讲Matlab求解微分方程(组)理论介绍：Matlab求解微分方程(组)命令求解实例：Matlab求解微分方程(组)实例实际应用问题通过数学建模所归纳得到的方程，绝大多数都是微分方程，真 正能得到代数方程的机会很少另一方面，能够求解的微分方程也是十分有限的， 特别是高阶方程和偏微分方程(组)这就要求我们必须研究微分方程(组)的 解法：解析解法和数值解法 一相关函数、命令及简介1.在Matlab中，用大写字母D表示导数，Dy表示y关于自变量的一阶导数， D2y表示y关于自变量的二阶导数，依此类推.函数dsolve用来解决常微分方程(组)的求解问题，调用

2、格式为：X=dsolve( eqn 1 , eqn2,)函数dsolve用来解符号常微分方程、方程组，如果没有初始条件，则求出通 解，如果有初始条件，则求出特解.注意，系统缺省的自变量为t2.函数dsolve求解的是常微分方程的精确解法，也称为常微分方程的符号解.但是，有大量的常微分方程虽然从理论上讲， 其解是存在的，但我们却无法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程数值解方面， MATLAB具有丰富的函数，我们将其统称为 solver，其一般格式为：T,Y=solver(odefu n, tspa n,yO)说明：(1)solver 为命令 ode45、ode23 od

3、e113ode15s ode23s ode23t、ode23tb、 ode15i 之一.(2)odefun是显示微分方程讨=f (t,y)在积分区间tspan = t,tf上从t0到tf用 初始条件y求解.(3)如果要获得微分方程问题在其他指定时间点 仏2,川,t f上的解，则令tspan=t。,t1,t2J|ltf(要求是单调的).(4)因为没有一种算法可以有效的解决所有的 ODE问题，为此，Matlab提供了多种求解器solver，对于不同的ODE问题，采用不同的solver.表1 Matlab中文本文件读写函数求解器ODE类型特点说明ode45非刚性单步算法：4、5阶Runge-Kutt

4、a方程；累计截断误差（Ax）3大部分场合的首选算法ode23非刚性单步算法：2、3阶Runge-Kutta方程；累计截断误差（Ax）3使用于精度较低的情形ode113非刚性多步法：Adams算法；高低精度可达10“10计算时间比ode45短ode23t适度刚性采用梯形算法适度刚性情形ode15s刚性多步法：Gears反向数值微分；精度中等若ode45失效时，可尝试使用ode23s刚性单步法：2阶Rosebrock算法；低精度当精度较低时，计算时间比ode15s短ode23tb刚性梯形算法；低精度当精度较低时，计算时间比ode15s短说明：ode23、ode45是极其常用的用来求解非刚性的标准形

5、式的一阶微分方程（组）的初值问题的解的Matlab常用程序，其中：ode23采用龙格-库塔2阶算法，用3阶公式作误差估计来调节步长，具有低 等的精度ode45则采用龙格-库塔4阶算法，用5阶公式作误差估计来调节步长，具有 中等的精度3.在matlab命令窗口、程序或函数中创建局部函数时，可用内联函数inline， inline函数形式相当于编写M函数文件，但不需编写 M-文件就可以描述出某种 数学关系 调用inline函数，只能由一个 matlab表达式组成，并且只能返回一个 变量，不允许u,v这种向量形式.因而，任何要求逻辑运算或乘法运算以求得最 终结果的场合，都不能应用inline函数，i

6、nline函数的一般形式为：FunctionName=inline（函数内容所有自变量列表例如：（求解F（x）=x2*cos（a*x）-b ,a,b是标量；x是向量）在命令窗口输入：Fofx=inline( x .A2*cos(a(*X) , x , a , b);g= Fofx(pi/3 pi/3.5,4,1)系统输出为：g=-1.5483 -1.7259注意：由于使用内联对象函数inline不需要另外建立m文件，所有使用比较 方便，另外在使用ode45函数的时候，定义函数往往需要编辑一个 m文件来单 独定义，这样不便于管理文件，这里可以使用 inline来定义函数.二实例介绍1.几个可以直

7、接用Matlab求微分方程精确解的实例例1求解微分方程y，2xy=xe程序：syms x y; y=dsolve( Dy+2*x*y=x*exp(-xA2)例2求微分方程xy y-ex =0在初始条件y=2e下的特解并画出解函数的图形.程序：syms x y; y=dsolve( x*Dy+y-exp(1)=0 y(1)=2*exp(1) ezplot(y)dx t5x y = e例3求解微分方程组(dt 在初始条件x|y = 1,y |t曲=0下的特解鱼-x3y=0.dt并画出解函数的图形.程序：syms x y tx,y=dsolve(,Dx+5*x+y=exp(t),Dy-x-3*y=0

8、,x(0)=1,y(0)=0,t,)simple(x);simple(y)ezplot(x,y,0,1.3);axis auto2.用ode23 ode45等求解非刚性标准形式的一阶微分方程(组)的初值问题的数值解(近似解)例4求解微分方程初值问题 dx = _2y 2 2x的数值解，求解范围为区【y(0)=i间0,0.5.程序：fun=inline(,-2*y+2*xA2+2*x,x,y,);x,y=ode23(fu n,0,0.5,1);plot(x,y,o-)例5求解微分方程马 7(1一丫2)巴 y =0, y(0) =1,y(0) =0 的解，并画出dt dt解的图形.分析：这是一个二

9、阶非线性方程，我们可以通过变换，将二阶方程化为一阶 方程组求解令Mg哼7，则捲(0) = 1x2(0) = 0dx!頁dX2 =7(1x2)X2 7, .dt编写M-文件vdp.mfun cti on fy=vdp(t,x) fy=x(2);7*(1-x(1)A2)*x (2) -x(1); end在Matlab命令窗口编写程序y0=1;0t,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折线

10、法求解的基本思想是将微分方程初值问题f (x, y).y(xO =y。化成一个代数(差分)方程，主要步骤是用差商y(x h) - y(x)替代微商 鱼，于是 h dx乂二 y(x。)记 Xk Xk h, yy(Xk),从而 yk y(Xk h),于是y。=y(x),/ x“=Xk+h, k = 0,1,2j|, n1yk 1 二 yk hf (Xk,yQ.例6用Euler折线法求解微分方程初值问题dy 2xy 2dx y2y(0)=1的数值解(步长h取0.4),求解范围为区间0,2.分析：本问题的差分方程为k =0,1,2,川，n-1=0, y0 = 1,h = 0.4xk 1 =Xk h,y

11、k 1 二 yk hf (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;% 数值解 for i=1: n-1y=y+h*subs(f,x,y,x,y);%subs，替换函数x=x+h;szj=szj;x,y;endszj 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

12、替换a和2替换b,返回cos(alpha)+sin(2)特别说明：本问题可进一步利用四阶 Runge-Kutta法求解，Euler折线法实际上就是一阶Runge-Kutta法，Runge-Kutta法的迭代公式为y。=y(x).Xk 厂Xk h.hyk 1 = yk (Li 2L2 2L3 L4),6k=0,1,2,|(, n-1Li 二 f (Xk, yk),h hL2 - f (Xk -,yk -Li),h hL f (Xk -,yk - L2),L4 = f(Xk h,yk hL3).相应的Matlab程序为： clear f=sym(y+2*X/yA2); a=0; b=2; h=0.

13、4; n=(b-a)/h+1; x=0; y=1; szj=x,y;% 数值解 for i=1: n-1l1= 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;endszj plot(szj(:,1),szj(:,2)练习与思考：(1)ode45求解问题并比较差异.利用Matlab求微分方程y-2y y = 0的解.H Q I(3)求解微

14、分方程 y 2(1-y)y y =0,0 乞 x 岂 30,y(0) hy(0) = 0 的特解., , Q H I I(4)利用Matlab求微分方程初值问题(1 x )y 2xy , y卜卫=1, y 1x=3的解.提醒：尽可能多的考虑解法 三微分方程转换为一阶显式微分方程组Matlab微分方程解算器只能求解标准形式的一阶显式微分方程(组)问题， 因此在使用ODE解算器之前，我们需要做的第一步，也是最重要的一步就是借 助状态变量将微分方程(组)化成 Matlab可接受的标准形式.当然，如果ODEs由一个或多个高阶微分方程给出，则我们应先将它变换成一阶显式常微分方程组 下面我们以两个高阶微分方程组构成的 ODEs为例介绍如何将它变换成一个一阶显式微分方程组.Step 1将微分方程的最高阶变量移到等式左边， 其它移到右边，并按阶次从 低到高排列.形式为：x(m)=f(t,x,x,x,川,x(m“y, y,yJ|,y(na)yn)=g(t,x,x,x,川xZyyyMyg)