微分方程数值格式的误差阶数计算新方法

刘小妹，陈涛

(中国传媒大学 理学院，北京 100024)

微分方程数值格式的误差阶数计算新方法

刘小妹，陈涛

(中国传媒大学 理学院，北京 100024)

针对常微分方程初值问题，在特定时间步长且真解未知的情况下，给出计算数值离散格式的误差随网格剖分尺寸的收敛阶数的新方法.通过数值实验，和必须已知真解情况下才能计算误差的传统方法进行对比分析.实验证明新方法可以不需要知道真解就能够对数值格式的误差阶数进行准确计算 突破了传统方法中的局限性.

常微分方程；数值求解；误差阶数

1 引言

在现代科学技术中，微分方程在很多领域中有着重要的应用，自动控制、浓度扩散、弹道的计算等。 常微分方程是描述连续变化的数学语言，在实际应用中，当时间或空间步长无限加密时，常微分方程的解收敛到原连续介质的真解。但对于特定步长时，在真解未知的情况下，人们很难得到常微分的精度做出可靠且实用的估计。针对这个问题，本文给出一种新方法，即真解未知的情况下分析常微分方程初值问题数值求解公式的阶数，根据实验真实可靠。

考虑一阶常微分方程的初值问题。

(1)

其中f为t和u的已知函数，u0为给定的初值。给予假设：设函数f(t，u)在区域：t0≤t≤T，|u|lt;∞内连续，并且关于u满足Lipschitz条件。可知，初值问题(1)在区间[t0，T]上有唯一解u(t)，并且u(t)是连续可微的。

2 数值解法

数值解法是一种离散化的方法，通过这种方法可以在一系列离散点t1，t2，…，tn上求出未知函数u(t)之值u(t1)，u(t2)，…，u(tN)的近似值u1，u2，…，uN。自变量t的离散值t1，t2，…，tN是事先取定的，通常取为等距的，即tn=t0+nh，n=1，2，…，N，其中hgt;0称为步长，而u1，u2，…，uN称为初值问题一个数值解。

为考察数值解法提供的数值解是否有实用价值，需要估计数值解与真解的误差，以便在实际计算中根据精度要求确定计算方案。之后还要知道，当步长充分小时，所得数值解能否足够精确地逼近初值问题的真解，即查看误差阶数，衡量计算格式随着网格尺寸较小的收敛速度。设un是在无舍入误差情况下，u(t)为真解，εn=‖u(tn)-un‖称为截断误差，若ε=Chp+Ο(hp+1)，则称该算法有p阶精度。

考虑初值问题(1)，本文给出几种常用的数值解法，其中包括显式Euler法、隐式Euler法、中点法、改进Euler法、梯形法等方法，并研究给出相应的截断误差的阶。

2.1 显式Euler法

显示Euler法的计算公式为

un+1=un+hf(tn，un).

此时截断误差

‖u(tn)-un‖=Ch+ο(h2).

其中C是系数，故显式Euler法的误差阶数为1，即当u(t)为一次多项式时，显式Euler是精确的。

2.2 隐式Euler法

隐式Euler法的计算公式为

此时同样截断误差

‖u(tn)-un‖=Ch+ο(h2)

故隐式Euler法为一阶方法。

2.3 中点法

中点法的计算公式为

un+1=un-1+2*hf(tn，un)

此时中点法的截断误差

‖u(tn)-un‖=Ch2+ο(h3)，

故中点法为二阶方法。

2.4 梯形法

梯形法的计算公式

此时同样截断误差

‖u(tn)-un‖=Ch2+ο(h3)，.

故梯形法为二阶方法。

2.5 改进Euler法

改进Euler法的计算公式

此时截断误差

‖u(tn)-un‖=Ch2+ο(h3)，.

故改进Euler法为二阶方法。

3 数值解估计误差新方法

新方法：设‖uk-u‖=Dτp+ο(τp+1)，其中Dgt;0是常数.可知‖uk-u‖区域零的速度与‖uk-uk+1‖的速度是相同的.那么则有下面公式

并且可得

(2)

在这种情况下我们可以消去真解u，所以如果利用这种方法，我们可以不用考虑真解就可以估计误差。

4 实例求解误差阶数

本文以一阶常微分方程为例子，利用MATLAB软件求解上述五种求数值解法的误差阶.我们分别以传统方法和我们的新方法来检验，并将两种结果做分析比较。

求解初值问题

这里我们取τ=1/25，则时间步长分别取1/25，1/26，1/27，1/28，1/29，1/210，1/211，1/212，1/213，1/214.

4.1 传统算法

图1 五种解法的回归结果

其得到的线性回归方程分别是：

y1=0.936490x+0.843358

y2=1.085501x+3.020986

y3=1.975451x+1.443158

y4=1.967441x+0.002437

y5=1.980239x+1.483173

其中yi和ri分别表示各种数值格式之间的误差离散点和相应误差回归直线。i=1为显示Euler法、i=2为隐式Euler法、i=3为中点法、i=4为梯形法、i=5为改进Euler法。

我们可以清晰地看出，在已知真解并给定时间步长的情况下，显示Euler法和隐式Euler法的误差阶为1，而中点法、梯形法、改进Euler法的误差阶为2，证明了结果和理论误差阶数吻合。

4.2 新方法

假设我们在不知道上述例子的真解是什么，传统方法无法做出误差分析，但按照本文给出求误差阶的新方法，我们得到五种求数值解方法的误差阶，结果如表1。

从表1我们看出，在未知真解并给定步长的情况下，显示Euler法、隐式Euler法的误差阶收敛为1，而中点法、梯形法、改进Euler法的误差阶收敛为2，和之前理论误差阶数吻合。另外，表2所展示的数据体现了误差阶p的收敛过程。实际计算误差时，我们只需要取最后三个网格尺寸，即k=8，9，10，相应的误差代入新方法的计算公式(2)，就可以很准确的计算误差阶数。

表1 五种解法的误差阶

以上两种方法做比较我们可以看出，两者得到的误差阶正如我们所期望的值是一样的，所以两种方法都真实可靠，但两者最大不同点是第一种已知真解，第二种真解未知。在实际应用中，大都很难得到真解，所以第二种，也就是本文提出的新方法在数值解法求误差精度时更具有应用性。

5 结论

本文在求解常微分初值问题求解误差阶时给出了一种新方法，即可以在真解未知的情况下可以得到数值解法的误差阶。利用五种数值解法，与已知真解的方法做出对比分析。分析结果表明：未知真解的情况下，通过验证和对比，本文给出的新方法可以得到数值解法的正确的精度阶，故新方法真实可靠。在实际应用中，新方法的最大优势在于不需要已知真解救能够对数值格式的误差阶数进行准确的估计，而传统方法必须计算多个网格剖分尺寸才能得到误差阶数。而且新方法不需要做回归，只需要计算三个网格尺寸的误差，即可得到误差阶数。

[1]胡建伟，汤怀民． 微分方程数值解法(第二版)[M]． 北京：科学出版社，2007．

[2]李庆扬，王能超，易大义． 数值分析(第五版)[M]． 北京：清华大学出版社，2008．

[3]戴嘉尊，邱建贤． 微分方程数值解法[M]． 南京：东南大学出版社，2002．

[4]谢焕田．应用Taylor公式分析常微分方程初值问题数值求解公式的精度[J].高师理科学刊，2011，(01)：9-11.

[5]魏明强. 一阶常微分方程数值解中四种算法的实例比较[J]．中国传媒大学学报(自然科学版)，2016，(04)：41-44.

(责任编辑：王 谦)

ANewMethodtoComputeErrorOrderofNumericalSchemeforDifferentialEquation

LIU Xiao-mei，CHEN Tao

(School of Science，Communication University of China，Beijing 100024，China)

A new method for computing the convergence order of the numerical scheme with respect to the size of the partition is given in the initial value problem of the first order ordinary differential equation and the exact solution is unknown. The new method is compared with the traditional method which should be implemented by known exact solution.The experimental results show that the new method can accurately obtain the error order of the numerical scheme even though the exact solution is unknown，and it can be used to estimate the error in the numerical scheme，and it is necessary to calculate the error in the numerical solution.

ordinary differential equation；numerical solution；error order

0241.1

A

1673-4793(2017)06-0044-05

2017-09-20

刘小妹(1993-)，女(汉族)，山东德州人, 中国传媒大学硕士研究生.E-mail：liuxiaomei@cuc.edu.com