常微分方程MATLAB解法
实验四 求微分方程的解
一、问题背景与实验目的
实际应用问题通过数学建模所归纳而得到的方程,绝大多数都是微分方程,真正能得到代数方程的机会很少.另一方面,能够求解的微分方程也是十分有限的,特别是高阶方程和偏微分方程(组).这就要求我们必须研究微分方程(组)的解法,既要研究微分方程(组)的解析解法(精确解),更要研究微分方程(组)的数值解法(近似解).
对微分方程(组)的解析解法(精确解),Matlab 有专门的函数可以用,本实验将作一定的介绍.
本实验将主要研究微分方程(组)的数值解法(近似解),重点介绍 Euler 折线法.
二、相关函数(命令)及简介
1.dsolve('equ1','equ2',„):Matlab 求微分方程的解析解.equ1、equ2、„为方程(或条件).写方程(或条件)时用 Dy 表示y 关于自变量的一阶导数,用用 D2y 表示 y 关于自变量的二阶导数,依此类推.
2.simplify(s):对表达式 s 使用 maple 的化简规则进行化简. 例如: syms x
simplify(sin(x)^2 + cos(x)^2) ans=1
3.[r,how]=simple(s):由于 Matlab 提供了多种化简规则,simple 命令就是对表达式 s 用各种规则进行化简,然后用 r 返回最简形式,how 返回形成这种形式所用的规则.
例如: syms x
[r,how]=simple(cos(x)^2-sin(x)^2) r = cos(2*x) how = combine
4.[T,Y] = solver(odefun,tspan,y0) 求微分方程的数值解. 说明:
(1) 其中的 solver为命令 ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb 之一.
dy
f(t,y)
(2) odefun 是显式常微分方程:dt
y(t0)y0
(3) 在积分区间 tspan=[t0,tf]上,从t0到tf,用初始条件y0求解.
(4) 要获得问题在其他指定时间点t0,t1,t2,上的解,则令 tspan= . [t0,t1,t2,,tf](要求是单调的)
(5) 因为没有一种算法可以有效地解决所有的 ODE 问题,为此,Matlab 提
供了多种求解器 Solver,对于不同的ODE 问题,采用不同的Solver.
(6) 要特别的是:ode23、ode45 是极其常用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的 Matlab 的常用程序,其中:
ode23 采用龙格-库塔2 阶算法,用3 阶公式作误差估计来调节步长,具有低等的精度.
ode45 则采用龙格-库塔4 阶算法,用5 阶公式作误差估计来调节步长,具有中等的精度.
5.ezplot(x,y,[tmin,tmax]):符号函数的作图命令.x,y 为关于参数t 的符号函数,[tmin,tmax] 为 t 的取值范围.
6.inline():建立一个内联函数.格式:inline('expr', 'var1', 'var2',…) ,注意括号里的表达式要加引号.
例:Q = dblquad(inline('y*sin(x)'), pi, 2*pi, 0, pi)
三、实验内容
1. 几个可以直接用 Matlab 求微分方程精确解的例子:
dyx2
2xyxe例1:求解微分方程,并加以验证. dx
求解本问题的Matlab 程序为:
syms x y %line1 y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') %line2 diff(y,x)+2*x*y-x*exp(-x^2) %line3 simplify(diff(y,x)+2*x*y-x*exp(-x^2)) %line4 说明:
(1) 行line1是用命令定义x,y为符号变量.这里可以不写,但为确保正确性,建议写上;
(2) 行line2是用命令求出的微分方程的解:
1/2*exp(-x^2)*x^2+exp(-x^2)*C1
(3) 行line3使用所求得的解.这里是将解代入原微分方程,结果应该为0,但这里给出:
-x^3*exp(-x^2)-2*x*exp(-x^2)*C1+2*x*(1/2*exp(-x^2)*x^2+exp(-x^2)*C1)
(4) 行line4 用 simplify() 函数对上式进行化简,结果为 0, 表明yy(x)的确是微分方程的解.
例2:求微分方程xy'yex0在初始条件y(1)2e下的特解,并画出解函数的图形.
求解本问题的 Matlab 程序为: syms x y
y=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x') ezplot(y)
eex
微分方程的特解为:y=1/x*exp(x)+1/x* exp (1) (Matlab格式),即y,
x
解函数的图形如图 1:
图1
dxt
5xyedt
例3:求微分方程组在初始条件x|t01,y|t00下的特解,
dyx3y0dt
并画出解函数的图形.
求解本问题的 Matlab 程序为: syms x y t
[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]);axis auto
微分方程的特解(式子特别长)以及解函数的图形均略. 2. 用ode23、ode45等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解(近似解).
dy2
2y2x2x
例4:求解微分方程初值问题dx的数值解,求解范围为区
y(0)1间[0, 0.5].
fun=inline('-2*y+2*x^2+2*x','x','y'); [x,y]=ode23(fun,[0,0.5],1); x'; y';
plot(x,y,'o-') >> x' ans =
0.0000 0.0400 0.0900 0.1400 0.1900 0.2400 0.2900 0.3400 0.3900 0.4400 0.4900 0.5000 >> y' ans =
1.0000 0.9247 0.8434 0.7754 0.7199 0.6764 0.6440 0.6222 0.6105 0.6084 0.6154 0.6179 图形结果为图 2.
图2
例 5:求解描述振荡器的经典的 Ver der Pol 微分方程
d2y2dy(1y)y0,y(0)1,y'(0)0,7.
dtdt2
分析:令x1y,x2
dx1dxdx2
,则1x2,2(1x1)x2x1. dtdtdt
先编写函数文件verderpol.m:
function xprime = verderpol(t,x) global mu;
xprime = [x(2);mu*(1-x(1)^2)*x(2)-x(1)]; 再编写命令文件vdp1.m: global mu; mu = 7; y0=[1;0]
[t,x] = ode45('verderpol',[0,40],y0); x1=x(:,1);x2=x(:,2); plot(t,x1)
图形结果为图3.
图3
3. 用 Euler 折线法求解
前面讲到过,能够求解的微分方程也是十分有限的.下面介绍用 Euler 折线法求微分方程的数值解(近似解)的方法.
Euler 折线法求解的基本思想是将微分方程初值问题
dy
f(x,y),
dx
y(x0)y0
化成一个代数方程,即差分方程,主要步骤是用差商于是:
y(xh)y(x)dy
替代微商,
hdx
记xk1
y(xkh)y(xk)
f(xk,y(xk)),
h
y0y(x0)
xkh,yky(xk),从而yk1y(xkh),则有
y0y(x0),
k0,1,2,,n1 xk1xkh,
yyhf(x,y).
kkkk1
例 6:用 Euler 折线法求解微分方程初值问题
2xdyy,2
y dx
y(0)1
的数值解(步长h取0.4),求解范围为区间[0,2].
解:本问题的差分方程为
x00,y01,h0.4,xk1xkh,
k0,1,2,,n1
2xyk1ykhf(xk,yk) (其中:f(x,y)y2).y
相应的Matlab 程序见附录 1. 数据结果为:
0 1.0000 0.4000 1.4000 0.8000 2.1233 1.2000 3.1145 1.6000 4.4593 2.0000 6.3074
图形结果见图4:
图4
特别说明:本问题可进一步利用四阶 Runge-Kutta 法求解,读者可将两个结果在一个图中显示,并和精确值比较,看看哪个更“精确”?(相应的 Matlab 程序参见附录 2).
四、自己动手
1. 求微分方程(x21)y'2xysinx0的通解. 2. 求微分方程y''2y'5yexsinx的通解. 3. 求微分方程组
dx
xy0dt
dyxy0dt
在初始条件x|t01,y|t00下的特解,并画出解函数yf(x)的图形. 4. 分别用 ode23、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解),求解区间为t[0,2].利用画图来比较两种求解器之间的差异.
5. 用 Euler 折线法求解微分方程初值问题
12x2y'y3,
y
y(0)1
的数值解(步长h取0.1),求解范围为区间[0,2].
6. 用四阶 Runge-Kutta 法求解微分方程初值问题
y'yexcosx,
y(0)1
的数值解(步长h取0.1),求解范围为区间[0,3].
四阶 Runge-Kutta 法的迭代公式为(Euler 折线法实为一阶 Runge-Kutta 法):
y0y(x0),
xk1xkh,h
yk1yk(L12L22L3L4)
6
k0,1,2,,n1 L1f(xk,yk)
hhL2f(xk,ykL1)
22
hhL3f(xk,ykL2)
22
L4f(xkh,ykhL3)
相应的 Matlab 程序参见附录 2.试用该方法求解第5题中的初值问题. 7. 用 ode45 方法求上述第 6 题的常微分方程初值问题的数值解(近似解),从而利用画图来比较两者间的差异.
五、附录
附录 1:(fulu1.m)
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]; for i=1:n-1
y=y+h*subs(f,{'x','y'},{x,y}); x=x+h;
szj=[szj;x,y]; end szj
plot(szj(:,1),szj(:,2))
附录 2:(fulu2.m)
clear
f=sym('y-exp(x)*cos(x)'); a=0; b=3; h=0.1;
n=(b-a)/h+1; x=0; y=1;
szj=[x,y]; for i=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))