常微分方程初值问题的MATLAB解法
大连圣亚海洋世界官网-
用
Matlab
求常微分方程
< br>(ODE)
的初值问题
(IVP)
本节考虑一阶常微分方程
u
f
(
t
,
u
< br>)
t
0
t
T
u
(
t
)
u
0
<
/p>
0
(1.1)
的数值求解问题，包括算
法公式及编程问题。对一阶常微分方程组的初值问题
u
1
(
t
0
)
u
1
0
f
1
(
t
,
u
1
,<
/p>
u
2
,
,
u
m
)
u
1
u
f
(
t
,
u
,
u
,
<
/p>
,
u
)
0
u
2
(
t
0
)
u
2
2
2
1
2
m
t
0
t
T
u
(
t
)<
/p>
u
0
f
m
(
t
,
u
1
,
u
2
,
,
u
m
)
u
m
m<
/p>
m
0
(1.2
)
可
以通过引入列向量
u
,
u
0
,
f
化成类似
(1.1)
的形式
其中
<
/p>
u
f
(
t
,
u
)
t
0
t
T
< br>
u
(
t
0
)
u
0
(1.
3)
u
1
(
t
)
<
/p>
u
1
(
t
0
)
f
1
(
t
,
u
1
,
u
2
,
,
u
m
)<
/p>
u
(
t
)
u
(
t
)
f
(
t
,
u
,
u
,
<
/p>
,
u
)
2
2
0
2
1
2
m
,
u
,
f
(
t
,
u<
/p>
)
u
(
t
)
0
<
/p>
u
(
t
)
u
(
t
)
f
(
t
,
u
,
u
,
,
u
)
1
2
m<
/p>
m
m
0
m
另外一个高阶常微分方程的初值问题
(1.4)
u<
/p>
(
m
)
f
(
t
,
u
,
u
,
u
,
,
u
(
m
1)
)
t
0
t
T
(
m
1)
(
m
1)
,
< br>
,
u
u
(
t
0
)
u
0
,
u
(
t
0
)
u
0
(
t
0
)
< br>
u
0
(1.5)
也可以通过变换
u<
/p>
1
u
,
u
2
u
,
u
3
u
,
,
u
m
u
(
m<
/p>
1)
化成维微分方程组：
< br>u
2
u
1
u
u
2
3
u
u
m
m
1
< br>f
(
t
,
u
1
,
u
2
,
,
u
m
)
u
m
(1.6)
我们
在设计算法时通常先对一维一阶常微分方程
(1.1)
进行，<
/p>
然后再将这个算法写成适合求解
(1.3)
的向量形式，
并以向量形式来进行编程。
1
、
非刚性系统与刚性系统
当
f
(
t
,
u
)
Au
(
t
)
时微分方程组
(1.3)
变成
如果系数矩阵<
/p>
A
R
m
m
u
Au
< br>(
t
)
特征值
1
,
2
,
,
m
满足
(1
.7)
1
2
m
，
对应的
特征向量为
v
1
,
v
2
,
,
v
m
,
则<
/p>
(1.7)
具
1
有通解
m
j
t
u
(
t
)
c
j
e
v
j
< br>
(
t
)
j
1
(
1.8)
其中
(
< br>t
)
为
(1.7)
的一个特解。如果
(1.7)
中的矩阵满足
(1) Re(
j
)
0
j
1,
2,
,
m<
/p>
(2)
s
Re(
1
)
Re(
m
)
1
(1.9)
就称微分方程组
(1.7)
是
刚性
的或坏条件的或病态的
。
s
称为
刚性比
。
对于一般
的一阶常微分方程组
(1.3)
，如果多元向量值函数
f
(
t
,
u
)
对
(1.4)
中向量
u
的分量的
Ja
cobi
矩阵
f
1
u
1
f
< br>2
J
(
t
)
u
1
f
m
u
1
f
< br>1
u
2
f
2
u
2
f
m
u
2
f
1
u
m
< br>
f
2
u
m
f
m
u
m
< br>
u
u
(
t
)
(1.10)
的特征值
1
(
t
),
2
(
< br>t
),
,
m
(
t
)
对于求解区间上的任何
t
[
t
0
,
T
]
满足
(1.9)
式，则称
微分方程组
(1.3)
为刚性的
。
2
、
解法器及调用格式
解法器
适用类型
Nonstiff
何时使用
使用算法
ode45
ode23
Most
of
the
time.
This
should
be
the
first
Runge-Kutta (4,5)
格式
solver you
try.
If
using
crude
error
tolerances
or
solving
Runge-Kutta (2,3)
格式
moderately
stiff problems.
If using
stringent error tolerances or solving
Adams-Bashforth-Moulton PECE
格式
a
computationally intensive ODE file.
If
ode45
is
slow
because
the
problem
is
stiff.
If using crude error tolerances to
solve stiff
systems and the
mass matrix is constant.
If
the problem is only moderately stiff and
you
need
a
solution
without
numerical
damping.
If using
crude error tolerances to solve stiff
systems.
Nonstiff
ode113
Nonstiff
ode15s
Stiff
ode23s
Stiff
ode23t
Moderately
Stiff
ode23tb
Stiff
有如下三种调用方法，其中
solver
是上
面三个解法器的名子。
(1) [t,y] =
solver(odefun,tspan,y0)
；
2
(2) [t,y] = so
lver(odefun,tspan,y0,options)
；
(3) [t,y] = solver(odefun,tspan,y0,o
ptions,p1,p2...)
；
它们输入参数是：
(1) odefun:
是一个函数句柄，指向初值问题
(1.3)
中的
函数
f
(
t
,
u
)
，它的基本形式是
dydt = f(t,y)
，
其中
dydt,y
是列
向量，而
t
是标量。
(2)
< br>tspan
：是一个向量，用于指定求解区间
[tspa
n(1),tspan(end)]
。要求这个向量的分量是有序的，或者单调上
升或者单调下降。
(3) y0
：是初值问题
(1.3)
中的初始列向量
u
0
。
(4) options
：用于对解法器缺省的求解参数进行设
置。如果不想改变缺省值，可以用空矩阵‘
［］
’来代替。
MATLAB
的
ODE
解法器设计了一个结构变量用于设置解法器的各项公共
参数，比如精度，
步长等。其中已
经定义了一些缺省值，用无参
数的
odeset
命令可以列出所有的公共参数及其缺省值：<
/p>
AbsTol: [ positive
scalar or vector {1e-6}
]
绝对误差
RelTol: [ positive scalar {1e-3}
]
相对误差
NormControl: [ on | {off}
]
是否用向量范数控制误差，缺省是按每个分量进行误差控制。
OutputFcn: [ function ]
OutputSel: [ vector of integers
]
Refine: [ positive integer
]
Stats: [ on | {off} ]
InitialStep: [ positive scalar
]
初始步长
MaxStep: [ positive scalar
]
最大步长，缺省为求解区间的十分之一。
BDF: [ on | {off} ]
是否利用向后差分来求
Jacobi
矩阵。
MaxOrder: [ 1 | 2 | 3 |
4 | {5} ]
Jacobian: [ matrix |
function ]
对刚性方程组的解法器可以提供相应的
J
acobi
矩阵信息
JPattern: [ sparse matrix ]
Vectorized: [ on | {off} ]
Mass: [ matrix | function ]
MStateDependence: [ none | weak |
strong ]
MvPattern: [ sparse
matrix ]
MassSingular: [ yes | no |
{maybe} ]
InitialSlope: [ vector ]
Events: [ function ]
每个解法器可以设置的参数如下表
:
Parameters
RelTol
,
AbsTol,
NormControl
OutputFcn
,
OutputSel
,
Refine
,
Stats
Events
MaxStep
,
InitialStep
Jacobian,
JPattern, Vectorized
ode45
ode23
ode113
ode15s
ode23s
ode23t
ode23tb
--
--
--
3
Mass
MStateDependence
MvPattern
MassSingular
InitialSlope
MaxOrder
,
BDF
--
--
--
--
--
--
--
--
--
--
--
--
--
--
--
--
--
--
--
--
--
(5)
输入参数
p1,p2...
等等是可选的参数
,
并且由解法器传递
给函数
odefun.
输出参数的含义是：
(1)
t
：节点列向量。
(2) y
：在节点处的解的近似值数组，第
i
行对应于
t
的第
i
行那
个节点处的近似解。即输出是如下形式：
节点
t
y
1
(t)
y
2
(t)
…………
y
m
(t)
例题
1<
/p>
、
分别用步长
h
0.05,0.1,0.125,0.2
来求初值问题在区间
[0,8]
上的近似解，
并与准确解<
/p>
u
1
2
t
进行图形上的比较。
t
u
u
<
/p>
2
u
0
t
8
u
(0)
1
(1.11)
解：利用如下代码，注意其中作为输入参数的函数的代码的编写方法，以及
ODE
解法器
ode45
的调用方法。
function zzz
h=0.1;
y0=1;tspan=0:h:8;
[t,y] =
ode45(@f,tspan,y0);
[t,y]
hold on;
plot(tspan,y,'--
k',tspan,fu(tspan),'-b');
title('Solution of IVP in eg1');
xlabel('time t');
ylabel('solution u');
legend('numerical solutons','true
solutions');
function
dudt=f(t,u)
dudt=u-2*t./u;
function u=fu(t)
u=sqrt(1+2*t);
由于结果较多我们只画出它们
的图形如图
1
。从图形中可以明显看出，当时间变量较大时比如
t
5
所求的
近似
解严重偏离了真解。
这说明初值问题
(1.11)
在时刻
0
有较好的稳定
性，
但是在时刻
4
以后可能稳定性就不
太好了，
为些我们利用
(1.11)
构
造一个初值问题：
t
u
< br>u
2
u
4
t
12
u
(4)
3
(1.12)
它与初值问题
(1.11)
< br>有相同的解。将上述代码适当改动来求解初值问题
(1.12)
< br>，并画出图形如图
2
，发现还是在初
始时刻
4
附近有较好的稳定性，所以我们在求解时不能希望
在很大的区间上求解。
4
图
1
图
2
例
2
、
求解二阶
Va
n der Pol
微分方程所构成的初值问题在区间
[0,
20]
上的近似解。
其中参数
是一个正数，
当较大时这个方程是刚性的。
u
(1
u
2
)
u
u
0 0
t
20
u
(0)
2,
u
(0)
0
u
< br>
v
2
v
(1
u<
/p>
)
u
u
0 0
t
20
u
(0)
2,
v
(0)
0
(1.13)
解：这是一个高阶方程，要用变量代换化成向量形式的一阶方程组形式
(1.3)
，令
v
u
则化成方程组
(1.14)
< br>T
T
2
v
,
(
1
u
)
v<
/p>
u
再令
u
(
t
)
u
(
t
),
v
(
t
)
，
f
(
t
,
u
)
<
/p>
则
(1.14)
的一阶向量形式为
(1.3)
。下列代码中编写了向量
值函数
f
(
t
,
u
)
，由还有一个可选参数
故要注意解法器的调用形
式。输出矩阵
u
的每一行代表解向量在每个节
< br>点处的近似函数值。对于本题来说输出矩阵
u
的第二列相
当于方程
(1.13)
解函数
u
(
t
)
导数的近似值
。可以对
0.5,1,2,5,1
0,100,1000
进行试验解法器
ode45
的求解结果及图像。
function zzz
%ode13.m
mu=10;mustr=num2str(mu);
h=0.1;
u0=[2,0]';tspan=0:h:20;
[t,u] = ode45(@f,tspan,u0,[],mu);
[t,u]
figure;
hold on;
plot(t,u(:,1),'--
k',t,u(:,2),'-b');
title(strcat('
例
2
中的初值问题
mu=',mustr));
xlabel('time
t');
ylabel('solution u');
l
egend('
数值解
','
导函数数
值解
');
hold off;
function dudt=f(t,u,mu)
5