第二讲 常微分方程数值求解

第二讲

常微分方程数值求解

—— MATLAB 求解

Matlab 解初值问题函数

 用 Maltab自带函数 解初值问题

 求解析解： dsolve

 求数值解：

ode45、ode23、

ode113、ode23t、ode15s、

ode23s、ode23tb

符号求解

dsolve

符号求解

dsolve 求解析解

 求解析解： dsolve

y=dsolve('eq1','eq2', ... ,'cond1','cond2', ... ,'v')

其中 y 为输出的解， eq1、eq2、... 为方程，cond1、

cond2、... 为初值条件，v 为自变量。

例 1： 求微分方程

的通解，并验证。

dx

+ = −

sol=dsolve('Dy+2*x*y=x*exp(-x^2)','x') syms x

diff(sol)+2*x*sol - x*exp(-x^2) % 验证

dsolve 的使用

 几点说明

 如果省略初值条件，则表示求通解；

 如果省略自变量，则默认自变量为

t

dsolve('Dy=2*x','x'); % dy/dx = 2x dsolve('Dy=2*x'); % dy/dt = 2x

 若找不到解析解，则提出警告，并返回空解。

 微分方程中用 D 表示对 自变量 的导数，如：

Dy y'； D2y y''； D3y y'''

dsolve 的使用

 使用符号方程

 导数：diff，如 diff(y)，diff(y,2)

 等号：==

 必须声明应变量与自变量!

例 1： 求微分方程

的通解，并验证。

dx

+ = −

syms y(x)

sol=dsolve(diff(y)+2*x*y==x*exp(-x^2)) diff(y)+2*x*y - x*exp(-x^2)

dsolve 举例

例 2：

下的特解，并画出解函数的图形。

0

'

xy + − y e =

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

ezplot(sol);

1 2 ( )

y = e

syms y(x)

sol=dsolve(diff(y)*x+y-exp(x)==0, y(1)==2*exp(1));

ezplot(sol);

syms y(x)

sol=dsolve(diff(y)*x+y-exp(x)==0, y(1)==2*exp(sym(1)));

ezplot(sol);

dsolve 举例

例3：

下的特解，并画出解函数的图形。

5

3 0 dx

x y e dt

dy x y

dt

 + + =

 

 − − =



[x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0', ...

'x(0)=1', 'y(0)=0', 't') ezplot(x,y,[0,1.3]);

0 0

1 0

|

|

t t

x y

=

=

 =

 =



注：解微分方程组时，如果所给的输出个数与方程个数相同，

则方程组的解按词典顺序输出；如果只给一个输出，则输出

dsolve 举例

例：

'x(0)=1', 'y(0)=1', 't')

sol = dsolve('Dx+5*x=0', 'Dy-3*y=0', ...

'x(0)=1', 'y(0)=1', 't')

这里返回的 sol 是一个 结构类型 的数据 sol.x % 查看解函数 x(t)

sol.y % 查看解函数 y(t)

dsolve的输出个数只能为一个或与方程个数相等

dsolve 举例

 使用符号方程

syms x(t) y(t)

sol=dsolve(diff(x)+5*x==0, diff(y)-3*y==0, ...

x(0)==1, y(0)==1)

dsolve 举例

dsolve('D2y=-a^2*y','y(0)=1','Dy(pi/a)=0')

例 4：

^{d y}

_{a y}

dt = −

0 1 π 0

( ) , '( / )

y = y a =

其中

a

syms y(t) a dy = diff(y);

sol=dsolve(diff(y,2)==-a^2*y, y(0)==1, dy(pi/a)==0)

数值求解

ode45、ode23、

ode113、ode23t、ode15s、

ode23s、ode23tb

数值求解

数值求解

[T, Y] = solver (odefun, tspan, y0)

其中 y0 为初值条件，tspan为求解区间；Matlab在数值求解 时自动对求解区间进行分割，T (列向量) 中返回的是分割点 的值(自变量)，Y (数组) 中返回的是这些分割点上的近似解，

其列数等于应变量的个数。

solver

ode113、ode15s、ode23s、ode23t、ode23tb）

没有一种算法可以有效地解决所有的 ODE 问题，因此

MATLAB 提供了多种ODE求解器，对于不同的ODE，可以 调用不同的求解器。

Matlab的ODE求解器

求解器 ODE类型 特点 说明

ode45 非刚性 单步法；4，5 阶 R-K 方法；

累计截断误差为 (△x)³

大部分场合的首选方法

ode23 非刚性 单步法；2，3 阶 R-K 方法；

累计截断误差为 (△x)³

使用于精度较低的情形

ode113 非刚性 多步法；Adams算法；高低精 度均可到 10^-3～10^-6

计算时间比 ode45 短

ode23t 适度刚性 采用梯形算法 适度刚性情形

ode15s 刚性 多步法；Gear’s 反向数值微 分；精度中等

若 ode45 失效时，可 尝试使用

ode23s 刚性 单步法；2 阶Rosebrock 算 法；低精度

当精度较低时，计算时 间比 ode15s 短

ode23tb 刚性 梯形算法；低精度 当精度较低时，计算时

ode15s

参数说明

odefun 为函数句柄，代表显式常微分方程，可以通过匿名

函数定义，或在函数文件中定义，然后通过函数句柄调用。

f = @(x,y) -2*y+2*x^2+2*x;

[x,y] = ode45(f, [0,0.5],1);

注：也可以在 tspan 中指定对求解区间的分割，如：

[x,y] = ode45(f, [0:0.1:0.5],1); % x=[0:0.1:0.5]

[T, Y] = solver (odefun,tspan,y0)

求 的数值解，求解区间为 [0,0.5]

2 2 2 2

(0) 1

dy y x x

dx y

 = − + +



 =



例 ：

数值求解举例

如果需求解的问题是高阶常微分方程，则需将其化为一阶常 微分方程组，此时必须用函数文件来定义该常微分方程组。

1

2

2 2

1 2 1

(1 )

7 dx x

dt

dx x x x

dt

µ µ

 =



 = − −



 = 令 ¹

2

x y x dy

dt

 =

 =



求解 Ver der Pol 初值问题

2

2

2 (1 ) 0

7, (0) 1, '(0) 0

d y dy

y y

dt dt

y y

µ µ

 − − + =



 = = =



例 ：

数值求解举例

 先编写函数文件 verderpol.m

function dx=verderpol(t,x) % 必须是两个输入和一个输出 global mu; % 其它参数只能通过全局变量传递

dx=[x(2); mu*(1-x(1)^2)*x(2) - x(1)]; % 必须是列向量

 然后编写脚本文件 vanderpol_main.m clear;

global mu;

mu=7;

y0=[1; 0];

[t,x] = ode45(@verderpol, [0,40], y0);

plot(t, x(:,1));