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

第二讲

常微分方程数值求解

—— 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 ：求微分方程 的通解，并验 证。sol=dsolve('Dy+2*x*y=x*exp(-x^2)','x')

syms x

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

2 ^x2

dy xy xe dx

  

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 ：求微分方程 的通解，并验 证。syms y(x)

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

2 ^x2

dy xy xe dx

  

dsolve

例 2 ：求微分方程 在初值条件 下的 特解，并画出解函数的图形。

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

ezplot(sol);

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);

0

' ^x

xy  y e  y( )1  2e

dsolve

例 3 ：求微分方程组 在初值条件

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

[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]);

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

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

5

3 0

dx t

x y e dt

dy x y

dt

   



   



0 0

1 0

|

|

t t

x y





 

 



dsolve

例：[x,y]=dsolve('Dx+5*x=0', 'Dy-3*y=0', ...

'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 ：求微分方程 的特解，初值条件为

其中

a

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

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

2

2 2

d y a y dt  

0 1  0

( ) , '( / )

y  y a 

数值求解

ode45 、 ode23 、

ode113 、 ode23t 、 ode15s

、

ode23s 、 ode23tb

数值求解

数值求解

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

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

solver

ode45 、 ode23 、 ode113 、 ode15s 、 ode23s 、 ode2 3t 、 ode23tb ）

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

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

Matlab

求解器 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 短 ode23t

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

间比 ode15s 短

参数说明

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

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

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

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

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

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

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

例 ： ^{2 2 2 2}

(0) 1

dy y x x

dx y

    



 



数值求解举例

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

令

求解 Ver der Pol 初值问题

例 ：

1

2

2 2

1 2 1

(1 )

7 dx x

dt

dx x x x

dt 



 



   



 

1

2

x y

x dy

dt

 

 



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));