matlab微分方程实例,Matlab常微分方程的解法

【实例简介】

和Matlab应用有关的，具体介绍常微分方程的使用和解法，原理性介绍，帮助理解。

局部截断误差指的是,按()式计算由到这一步的计算值与精确值

之差

+。为了估计它,由

展开得到的精确值

是

()、()两式相减(注意到=

)得

即局部截断误差是阶的,而数值算法的精度定义为:

若一种算法的局部截断误差为

,则称该算法具有阶精度。

显然葳大,方法的精度葳高。式()说明,向前

方法是一阶方法,因此

它的精度不高。

改进的方法

梯形公式

利用数值积分方法将微分方程离散化时,若用梯形公式计算式中之右端积分,

即

并用代替

,则得计算公式

这就是求解初值问题()的梯形公式。

直观上容易看出,用梯形公式计算数值积分要比矩形公式好。梯形公式为二阶方法

梯形公式也是隐式格式,一般需用迭代法求解,迭代公式为

由于函数

关于满足

条件,容易看出

其中为

常数。因此,当

如果实际计算时精度要求不太高,用公式()求解时,每步可以只迭代一次,由此导

出一种新的方法一改进

改进

法

按式()计算问题()的数值解时,如果每步只迭代一次,相当于将公式

与梯形公式结合使用:先用公式求的个初步近似值,称为预测值,然

后用梯形公式校正求得近似值,即

预测

校正

式()称为由公式和梯形公式得到的预测一校正系统,也叫改进法。

为便于编制程序上机,式()常改写成

改进法是二阶方法

龙收一库塔(

-)方法

回到方法的基本思想—用差商代替导数一上来。实际上,按照微分中佰定理

应有

+b<6<

注意到方程=

就有

不妨记

+日,称为区问

上的平均斜率。可见给出一种

斜率,()式就对应地导出一种算法

向前公式简单地取

为,精度自然很低。改进的公式可理

解为取

的平均值,其中

,这种处

理提高了精度

如上分析启示我们,在区间

内多取几个点,将它们的斜率加权平均作为

就有可能构造出精度史高的计算公式。这就是龙格一库塔方法的基本思想。

首先不妨在区间内仍取个点,仿照()式用以下形式试一下

+2+a

+a+B

其中元λαβ为待定系数,看看如何确定它们使()式的精度尽量高。为此我们

分析局部截断误差

因为

,所以()可以化为

182

2+2

+a

B

其中在点

作了

展开。()式又可表为

+元+2

注意到

+,可见为使误差

只须令

+元

c

待定系数满足()的()式称为阶龙格一库搭公式。由于()式有个未知数

而只有个方程,所以解不唯一。不难发现,若令A=元

a=B=,即为改

进的

公式。可以证明,在

内只取点的龙格一库塔公式精度最高为

阶

阶龙格一库塔公式

要进一步提高精度,必须取更多的点,如取点构造如下形式的公式:

+++1+

+a

+B +B+B

其中待定系数λαβ共个,经过与推导阶龙格一库塔公式类似、但更复杂的计

算,得到使局部误差

的个方程。取既满足这些方程、又较简

单的组2aB,可得

这貮是常用的阶龙格一库塔方法(简称方法)。

线性多步法

以上所介绍的各种数值解法都是单步法,这是因为它们在计算时,都只用到

前步的值,单步法的般形式是

其中φ

称为增量函数,例如方法的增量函数为

,改进法的

增量函数为

如何通过较多地利用前面的已知信息,如

,来构造高精度的算法

计算,这就是多步法的基本思想。经常使用的是线性多步法

让我们试验一下

即利用

计算的情况。

从用数值积分方法离散化方稈的()式

出发,记

,式中被积函数

用二节点

的插值公式得到(因

,所以是外插

此式在区间

上积分可得

于是得到

注意到插值公式()的误差项含因子

,在区间

上积分后

181

得出,故公式()的局部截断误差为

,精度比向前公式提高阶。

若取=

可以用类似的方法推导公式,如对于=有

其局部截断误差为

如果将上面代替被积函数

用的插值公式由外插改为内插,可进一步减

小误差。内插法用的是

取三吋得到的是梯形公式,取=时

可得

与()式相比,虽然其局部截断误差仍为

,但因它的各项系数(绝对值)大

为减小,误差还是小了。当然,()式右端的未知,需要如同向后公式一

样,用迭代或校正的办法处理。

阶微分方程组与高阶微分方程的数值解法

阶微分方程红的数值解法

设有一阶微分方程组的初值问题

若记

,则初值

问题()可写成如下冋量形式

如果向量函数

在区域

连续,且关于满足

条件,即存在

使得对∈

都有

那么问题()在上存在唯解

问题()与()形式上完全相同,故对初值闩题()所建立的各种数值解法可

仝部用于求解问题()。

高阶微分方程的数值解法

高阶微分方稈的初值问题可以通过变量代换化为一阶徼分方程组初值问题

设有阶常微分方程初值问题

引入新变量=

,问题()就化为一阶微分方程初值问题

然后用中的数值方法求解问题(),就可以得到问题()的数值解。

最后需要指岀的是,在化学工程及自动控制等领域中,所涉及的常微分方程组初值

问题常常是所谓的“刚性”问题。具体地说,对·阶线性微分方程组

其中Φ∈

为阶方阵。若知阵的特征值λ

满足关系

>>

<

<<

则称方程组()为刚性方程组或方程组,称数

<<

≤

为刚性比。对刚性方程组,用前面所介绍的方法求解,都会遇到木质上的困难,这是由

数值方法本身的稳定性限制所决定的。理论上的分析表明,求解刚性问题所选用的数值

方法最好是对步长不作任何限制。

解法

数值解

非刚性常微分方程的解法

的工具箱提供了几个解非刚性常微分方程的功能函数,如,

,其中采用四五阶方法,是解非刚性常微分方程的首选方法,

采用二三阶方法,

采用的是多步法,效率一般比

同

的工具箱中没有

方法的功能函数。

()对简单的一阶方稈的初值问题

改进的公式为

我们自己编写改进的方法函数

如下:

186

例用改进的方法求解

解编写函数文件

如下

在 Matlab命令窗口输入

x, y]=eulerpro( doty', 0, 0.5, 1, 10)

即可求得数值解。

(TT)ode23,ode45,ode13的使用

Matlab的函数形式是

J=solver(F', tspan, yO

这里ver为ode45,ode23,ode113,输入参数F是用M文件定义的微分方程

右端的函数。 tspan-[t0, tfinal]是求解区间,y0是初值。

例2用RK方法求解

解同样地编写函数文件

如下

在Mat1ab命令窗口输入:

x,y]-ode45(doty3,0,0.5,1)

即可求得数值解。

刚性常微分方程的解法

的工具箱提供了几个解刚性常微分方程的功能函数,如

,这些函数的使用同上述非刚性微分方程的功能函数。

高阶微分方程

解法

例考虑初值问题

解()如果设

,那么

初值问题可以写成=

=的形式,其中

()把一阶方程组写成接受两个参数和,返回一个列向量的文件

注意:尽管不一定用到参数和

文件必须接受此两参数。这里向量必须是列

向量。

(ii)用 Matlab解决此问题的函数形式为

[T, Y=salver('F, tspan, yo)

这里 solver为ode45、odc23、odel13,输入参数F是用M文件定义的常微分方程组,

tspan-[ to trina]是求解区间,y0是初值列向量。在 Matlab命令窗口输入

[T,Y]=ode45(F,[01],[0;1;1])

扰得到上述常微分方程的数值解。这里Y和时刻T是一一对应的,Y(:,1)是初值问题的

解,Y(∷,2)是解的导数,Y(:,3)是解的二阶导数

例求 van der pol方程

的数值解,这里>是·参数。

解(i)化成常微分方程组。设

则有

()书写文件(对于=)

()调用

函数。对于初值

=,解为

()观察结果。利用图形输出解的结果:

例 van der pol方程,=

(刚性)

鲆(i)书写M文件vdp10.m

188

【实例截图】

【核心代码】