常微分方程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 求微分方程精确解的例子：

dyx2

2xyxe例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， 表明yy(x)的确是微分方程的解．

例2：求微分方程xy'yex0在初始条件y(1)2e下的特解，并画出解函数的图形．

求解本问题的 Matlab 程序为： syms x y

y=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x') ezplot(y)

eex

微分方程的特解为：y=1/x*exp(x)+1/x* exp (1) (Matlab格式)，即y，

x

解函数的图形如图 1：

图1

dxt

5xyedt

例3：求微分方程组在初始条件x|t01,y|t00下的特解，

dyx3y0dt

并画出解函数的图形．

求解本问题的 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

2y2x2x

例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(1y)y0,y(0)1,y'(0)0,7.

dtdt2

分析：令x1y,x2

dx1dxdx2

,则1x2,2(1x1)x2x1. 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(xh)y(x)dy

替代微商,

hdx

记xk1

y(xkh)y(xk)

f(xk,y(xk)),

h

y0y(x0)

xkh,yky(xk)，从而yk1y(xkh)，则有

y0y(x0),

k0,1,2,,n1 xk1xkh,

yyhf(x,y).

kkkk1

例 6：用 Euler 折线法求解微分方程初值问题

2xdyy,2

y dx

y(0)1

的数值解(步长h取0.4)，求解范围为区间[0,2]．

解：本问题的差分方程为

x00,y01,h0.4,xk1xkh,

k0,1,2,,n1 

2xyk1ykhf(xk,yk) （其中：f(x,y)y2).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. 求微分方程(x21)y'2xysinx0的通解． 2. 求微分方程y''2y'5yexsinx的通解． 3. 求微分方程组

dx

xy0dt



dyxy0dt

在初始条件x|t01,y|t00下的特解，并画出解函数yf(x)的图形． 4. 分别用 ode23、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解)，求解区间为t[0,2]．利用画图来比较两种求解器之间的差异．

5. 用 Euler 折线法求解微分方程初值问题

12x2y'y3,

y

y(0)1

的数值解(步长h取0.1)，求解范围为区间[0,2]．

6. 用四阶 Runge-Kutta 法求解微分方程初值问题

y'yexcosx,



y(0)1

的数值解(步长h取0.1)，求解范围为区间[0,3]．

四阶 Runge-Kutta 法的迭代公式为(Euler 折线法实为一阶 Runge-Kutta 法)：



y0y(x0),

xk1xkh,h

yk1yk(L12L22L3L4)

6

k0,1,2,,n1 L1f(xk,yk)

hhL2f(xk,ykL1)

22

hhL3f(xk,ykL2)

22

L4f(xkh,ykhL3)

相应的 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))

实验四 求微分方程的解

一、问题背景与实验目的

实际应用问题通过数学建模所归纳而得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少．另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）．这就要求我们必须研究微分方程（组）的解法，既要研究微分方程（组）的解析解法（精确解），更要研究微分方程（组）的数值解法（近似解）．

对微分方程（组）的解析解法(精确解)，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 求微分方程精确解的例子：

dyx2

2xyxe例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， 表明yy(x)的确是微分方程的解．

例2：求微分方程xy'yex0在初始条件y(1)2e下的特解，并画出解函数的图形．

求解本问题的 Matlab 程序为： syms x y

y=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x') ezplot(y)

eex

微分方程的特解为：y=1/x*exp(x)+1/x* exp (1) (Matlab格式)，即y，

x

解函数的图形如图 1：

图1

dxt

5xyedt

例3：求微分方程组在初始条件x|t01,y|t00下的特解，

dyx3y0dt

并画出解函数的图形．

求解本问题的 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

2y2x2x

例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(1y)y0,y(0)1,y'(0)0,7.

dtdt2

分析：令x1y,x2

dx1dxdx2

,则1x2,2(1x1)x2x1. 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(xh)y(x)dy

替代微商,

hdx

记xk1

y(xkh)y(xk)

f(xk,y(xk)),

h

y0y(x0)

xkh,yky(xk)，从而yk1y(xkh)，则有

y0y(x0),

k0,1,2,,n1 xk1xkh,

yyhf(x,y).

kkkk1

例 6：用 Euler 折线法求解微分方程初值问题

2xdyy,2

y dx

y(0)1

的数值解(步长h取0.4)，求解范围为区间[0,2]．

解：本问题的差分方程为

x00,y01,h0.4,xk1xkh,

k0,1,2,,n1 

2xyk1ykhf(xk,yk) （其中：f(x,y)y2).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. 求微分方程(x21)y'2xysinx0的通解． 2. 求微分方程y''2y'5yexsinx的通解． 3. 求微分方程组

dx

xy0dt



dyxy0dt

在初始条件x|t01,y|t00下的特解，并画出解函数yf(x)的图形． 4. 分别用 ode23、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解)，求解区间为t[0,2]．利用画图来比较两种求解器之间的差异．

5. 用 Euler 折线法求解微分方程初值问题

12x2y'y3,

y

y(0)1

的数值解(步长h取0.1)，求解范围为区间[0,2]．

6. 用四阶 Runge-Kutta 法求解微分方程初值问题

y'yexcosx,



y(0)1

的数值解(步长h取0.1)，求解范围为区间[0,3]．

四阶 Runge-Kutta 法的迭代公式为(Euler 折线法实为一阶 Runge-Kutta 法)：



y0y(x0),

xk1xkh,h

yk1yk(L12L22L3L4)

6

k0,1,2,,n1 L1f(xk,yk)

hhL2f(xk,ykL1)

22

hhL3f(xk,ykL2)

22

L4f(xkh,ykhL3)

相应的 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))