微分方程数值解法李荣华答案

微分方程数值解法李荣华答案

【篇一：高阶常微分方程的数值求解】

t>谷照升

（长春工程学院理学院，长春，130012）

摘要对经典初始条件的高阶常微分方程，给出其数值求解方法。该方法比runge-kutta法具有更好的适应性、易用性、计算速度和可控制的更高精度。

关键词常微分方程；数值解；算法

中图分类号: o241.81文献标识码: a 1

1 引言

求解复杂的1阶常微分方程，通常只能采用数值解法。数值解法一般又以runge-kutta法为主。对高阶常微分方程，则通常是将其转化为1阶常微分方程组，再用runge-kutta法求解[1, 2]。但这种通用性方法，在精度、软件计算的适应度方面，常常不够理想，甚至得不到结果。针对不同的广普性方程类型，可以建立更具针对性的计算方法。例如，针对三类边值条件和特定形式的方程，已经有有相应的差分法和有限元法。每一种更具针对性的方法，都有其更高的精度和更为健壮的算法，当然也存在其必然的局限性。

本文对如下形式的2阶常微分方程

?d2ydy?a(x)?b(x)y?f(x)?2dxdx?? x∈[a, b]

(1) ?y(a)?ya?y?(a)?y?a???

给出了完整、方便、高效的单步法数值求解算法，并将其推广到任意高阶问题的求解。 2 算法

思路：将[a, b]分割为a=x0 x1…xn=b，步长?xi?xi?xi?1。从i=1开始，利用数值积分公式，通过[xi-1, xi]上的数值积分，先求

y??(xi)，再求出y?(xi)，最后求y(xi)。依次取i=2, 3, … ,n，得到各点y??(xi)、y?(xi)、y(xi)的近似值。

根据数值积分的算法不同，主要有两种不同的计算公式。

2.1、矩形数值积分算法

矩形算法形式简单，精度偏低。基于不同的积分形式，又可以得到3种不同的计算公式，其精1 基金项目：吉林省自然科学基金项目（201215115）

1

度无显著差别。此处仅给出与2.2梯形算法最为接近的一种。

利用y?(xi)?y?(xi?1)??xi

xi?1y??(x)dx?y??(xi)?xi

y(xi)?y(xi?1)??

在xi点得到 xixi?1y?(x)dx?y?(xi)?xi

y?(xi)?y??(xi)?xi?y?(xi?1) (2)

2y(xi)?y?(xi)?xi?y(xi?1)?y??(xi)?xi?y?(xi?1)?xi?y(xi?1)(3)

这里 i?1,2,...,n。

将(2) 、(3)代入(1), 在xi点得到

2y??(xi)?a(xi)[y??(xi)?xi?y?(xi?1)]?b(xi)[y??(xi)?xi?y?(xi?1)?xi ?y(xi?1)]?f(xi) ???xi?y(xi?1)] (4) 从而

y??(xi)?f(xi)?a(xi)y(xi?1)?b(xi)[y(xi?1)

1?a(xi)?xi?b(xi)?xi

算法：

step 1: i=1;

?代入 (4), 输出y??(x1); 将(1)中 y(x0)?ya、y?(x0)?ya

将y??(x1)代入 (2), 输出y?(x1);

取y(x1)?y?(x0)?x1?y(x0), 输出y(x1)。

step 2: while in

{

i= i+1;

将y(xi?1)、y?(xi?1)代入 (4), 输出y??(xi);

将y??(xi)代入 (2), 输出y?(xi);

取y(xi)?y?(xi?1)?xi?y(xi?1), 输出y(xi);

}

2.2、梯形数值积分算法

根据(1)式，首先算出

y??(x0)?f(x0)?a(x0)y?(x0)?b(x0)y(x0)(5)

2

对i?1,2,3,...,n, 取y?(xxi

i)?y?(xi?1)??xy??(x)dx?1[y??(xi)?y??(xi?1)]?xi,

i?12

y(xxi1

i)?y(xi?1)??xy?(x)dx?[y?(xi)?y?(x

?12i?1)]?xi,

i

则xi点满足

y?(x1

i)?2[y??(xi)?y??(xi?1)]?xi?y?(xi?1) (6)

y(x1

i)?2[y?(xi)?y?(xi?1)]?xi?y(xi?1)

?1

2{1

2[y??(xi)?y??(xi?1)]?xi?y?(xi?1)?y?(xi?1)}?xi?y(xi?1) (7)

?1

4y??(x21

4??(x2

i)?xi?yi?1)?xi?y?(xi?1)?xi?y(xi?1),

将(6) 、(7)代入(1), 得到

y??(x?a(x1

i)i2[y??(xi)?y??(xi?1)]?xi?y?(xi?1)}

?b(x121

2

i)[4y??(xi)?xi?4y??(xi?1)?xi?y?(xi?1)?xi?y(xi?1)]?f(xi)

从而

f(xa(x1(x12

y??(xi)?i)[y??(xi?1)?xi?y?i?1)]?b(xi)[y??(xi?1)?xi?y?(xi?1)?xi? y(x

i)?24i?1)]

1?1

2a(x12

i)?xi?4b(xi)?xi

i?1,2,3,...,n

算法：

step 1:

按(5)式求出y??(x0)。

step 2:

i=0; while in

{

i= i+1;

将y(xi?1)、y?(xi?1)、y??(xi?1)代入 (8), 输出y??(xi);

3 (8)

将y??(xi?1)、y??(xi)代入 (6), 输出y?(xi);

将y(xi?1)、y?(xi?1)、y?(xi)代入 y(xi)?

} 1[y?(xi)?y?(xi?1)]?xi?y(xi?1), 输出y(xi); 2

3 误差分析