微分方程數值格式的誤差階數計算新方法

劉小妹，陳濤

(中國傳媒大學 理學院，北京 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