大连圣亚海洋世界官网-2021年2月7日发(作者：转身之后还是你)

用

Matlab

求常微分方程

<

br>(ODE)

的初值问题

(IVP)

本节考虑一阶常微分方程



u





f

(

t

,

u

<

br>)

t

0



t



T



u

(

t

)



u

0

<

/p>

0

(1.1)

的数值求解问题，包括算

法公式及编程问题。对一阶常微分方程组的初值问题



u

1

(<

p>

t

0

)



u

1

0





f

1

(

t

,

u

1

,<

/p>

u

2

,



,

u

m

)<

p>



u

1





u





f

(

t

,

u

,

u

,

<

/p>

,

u

)

0

u

2

(

t<

p>

0

)



u

2



2



2

1

2

m

t

0



t



T

















u

(

t

)<

/p>



u

0





f

m

(<

p>

t

,

u

1

,

u

2

,



,

u

m

)



u

m

m<

/p>



m

0

(1.2

)







可

以通过引入列向量

u

,

u

0

,

f

化成类似<

p>

(1.1)

的形式

其中



<

/p>







u





f

(<

p>

t

,

u

)

t

0



t



T



<

br>







u

(

t

0

)



u

0

(1.

3)



u

1

(

t

)

<

/p>



u

1

(

t

0

)

<

p>



f

1

(

t

,

u

1

,

u

2

,



,

u

m

)<

/p>















<

p>

u

(

t

)

u

(

t

)

f

(

t

,

u

,

u

,

<

/p>

,

u

)







2

2<

p>

0

2

1

2

m



,

u







,

f

(

t

,

u<

/p>

)







u

(

t<

p>

)





0

























<

/p>





u

(

t

)

u

(<

p>

t

)

f

(

t

,

u

,

u

,



,

u

)

1

2

m<

/p>





m





m

0

<

p>



m

另外一个高阶常微分方程的初值问题

(1.4)

u<

/p>

(

m

)



f

(

t

,<

p>

u

,

u



,

u





,



,

u

(

m



1)

)

t

0



t

p>



T





(

m



1)

(

m



1)





,

<

br>

,

u

u

(

t

0

)



u

0

,

u

p>



(

t

0

)



u

0

(

t

0

)

<

br>

u

0





(1.5)

也可以通过变换

u<

/p>

1



u

,

u

2



u<

p>



,

u

3



u





,



,

u

m



u

(

m<

/p>



1)

化成维微分方程组：





<

br>u

2



u

1



u





u

2

3



p>







p>



u





u

m



m



1





<

br>f

(

t

,

u

1

,

u

2

,



,

u

p>

m

)





u

m

(1.6)

我们

在设计算法时通常先对一维一阶常微分方程

(1.1)

进行，<

/p>

然后再将这个算法写成适合求解

(1.3)

的向量形式，

并以向量形式来进行编程。

1

、

非刚性系统与刚性系统

<

p>









当

f

(

t

,

u

)



Au





(

t

)

时微分方程组

(1.3)<

p>

变成

如果系数矩阵<

/p>

A



R

m



m

<

p>





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



p>

u

(

t

)





c

j

e

v

j



<

br>

(

t

)

j



1

(

1.8)

其中



(

<

br>t

)

为

(1.7)

的一个特解。如果

(1.7)

中的矩阵满足



(1) Re(



j

)



p>

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

矩阵







p>





f

1





u



1





f

<

br>2



J

(

t

)







u

1





p>







f

m





u



1



f

<

br>1



u

2



f

2



u

2





f

p>

m



u

2



f

1





u

m



<

br>



f

2







u

m







p>









f

m







u

m



<

br>





u



u

(

t

)



(1.10)

的特征值<

p>



1

(

t

),



2

(

<

br>t

),



,



m

(

t

)

对于求解区间上的任何

t



[

t

0

,

T

]

满足

(1.9)

式，则称

微分方程组

(1.3)

为刚性的

。

2<

p>

、

解法器及调用格式

解法器

适用类型

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<

p>

：是初值问题

(1.3)

中的初始列向量

u

0

。

(4) options

：用于对解法器缺省的求解参数进行设

置。如果不想改变缺省值，可以用空矩阵‘

［］

’来代替。

p>

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

矩阵。<

p>

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

p>





u

(0)

p>



1

(1.11)

解：利用如下代码，注意其中作为输入参数的函数的代码的编写方法，以及

ODE

p>

解法器

ode45

的调用方法。<

p>

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

例

p>

2

、

求解二阶

Va

n der Pol

微分方程所构成的初值问题在区间

[0,

20]

上的近似解。

其中参数



是一个正数，

当较大时这个方程是刚性的。



u



p>







(1



u

2

)<

p>

u





u



0 0



t



20





u

(0)



2,

u



(0)



0



u

<

br>



v



2



v







(1



u<

/p>

)

u





u



0 0



p>

t



20



u

(0)



2,

v

(0)



p>

0



(1.13)

解：这是一个高阶方程，要用变量代换化成向量形式的一阶方程组形式

(1.3)

p>

，令

v



u



则化成方程组

(1.14)





<

br>T



T

2



v

,



(

1



u

)

v<

/p>



u

再令

u

p>

(

t

)





u

(

t

),

v

(

t

)



，

f

(

t

,

u

)







<

/p>

则

(1.14)

的一阶向量形式为

(1.3)

。下列代码中编写了向量





值函数

f

(

t

,

u

)<

p>

，由还有一个可选参数



故要注意解法器的调用形

式。输出矩阵

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

大连圣亚海洋世界官网-

大连圣亚海洋世界官网-

大连圣亚海洋世界官网-

大连圣亚海洋世界官网-

大连圣亚海洋世界官网-

大连圣亚海洋世界官网-

大连圣亚海洋世界官网-

大连圣亚海洋世界官网-