基于MATLAB一階微分方程的仿真

任春年 李會榮 魏倩茹

摘? 要：為得到更加精確的求解微分方程數(shù)值解的方法，采用了顯式歐拉法、隱式歐拉法、改進歐拉法以及四階龍格庫塔等方法與MATLAB軟件中專有的ode45函數(shù)作比較，對用不同方法來求解微分方程的求解結(jié)果進行了研究，通過例證以及數(shù)據(jù)分析，得出在步長h任意時，四階龍格庫塔法的精準(zhǔn)度、穩(wěn)定性都要高于其他三種歐拉法，使在微分方程求解方法的選擇上更具針對性。

關(guān)鍵詞：微分方程;歐拉法;四階龍格庫塔法

中圖分類號：TP391.9? ? ? ? ? ? ? ?文獻標(biāo)識碼：A文章編號：2096-4706（2021）15-0160-03

Abstract： In order to obtain a more accurate numerical solution of differential equations， the explicit Euler method， implicit Euler method， improved Euler method and fourth-order Runge-Kutta method are compared with the special ode45 function of MATLAB software. The results of solving differential equations by different methods are studied. Through examples and data analysis， when the step h is arbitrary， the accuracy and stability of the fourth-order Runge-Kutta method are higher than the other three Euler methods， which makes the choice of differential equation solution method more targeted.

Keywords： differential equation; Euler method; fourth-order Runge-Kutta method

0? 引? 言

微分方程的理論和方法不僅廣泛應(yīng)用于自然科學(xué)，在數(shù)學(xué)研究領(lǐng)域內(nèi)，通過深研相關(guān)文獻可知，關(guān)于一階微分方程的各種邊值問題已取得不少的研究成果。文獻[1-2]中表示歐拉法解微分方程時過程雖簡單，但是誤差較大。

本文利用顯式歐拉法、隱式歐拉法、改進歐拉法、四階龍格庫塔等方法來對一階微分方程進行求解，通過對比相關(guān)數(shù)據(jù)，分析以上各種方法求解一階微分方程的優(yōu)劣勢，進一步提高一階微分方程的求解效率。

1? 常微分方程的基本概念

含有未知函數(shù)導(dǎo)數(shù)（或微分）的方程，稱為微分方程，未知函數(shù)是一元函數(shù)的微分方程稱為常微分方程，未知函數(shù)是多元函數(shù)的稱為偏微分方程[3]。

一階常微分方程初始問題：

其中x0為初始時間（已知常數(shù)），y0為初始狀態(tài)，函數(shù)f（x，y）關(guān)于y滿足利普希茨條件，保證初值問題解的存在且唯一。

2? 歐拉法

歐拉法的主要思想是迭代思想，也就是逐次替代，從而求出所需要的解，而且還要達到一定的精準(zhǔn)度。歐拉法分別有顯式歐拉法（前向歐拉法）、隱式歐拉法（后向歐拉法）和改進歐拉法。如果微分方程的最本質(zhì)的特征是方程中含有倒數(shù)項，那么數(shù)值解首要先消除倒數(shù)值，這個過程被稱為離散化。

2.1? 顯式歐拉法

使用一階前向差商代替微分[4]：

式中，h為時間步長。

微分方程的顯式差分方程：

上式是關(guān)于y（m）向y（m+1）的遞推形式，則可逐步求出y（1），y（2），y（3）…y（n）…，是關(guān)于y（m）的一個直接計算公式，既所求出的離散序列為該方程的數(shù)值解。

2.2? 隱式歐拉法

使用一階向后差商代替微分[5]：

微分方程的隱式差分方程：

在第m步迭代時，y（m-1），x（m）已知，y（m）為待求未知量，f（x（m），y（m））是關(guān)于待求未知量y（m）的函數(shù)，則方程是關(guān)于y（m）的非線性方程。通過非線性方程的迭代解法求出y（m）。則可逐步求出y（1），y（2），y（3）…y（n）…，所求出的離散序列為方程的數(shù)值解。

2.3? 改進歐拉法

先用顯式歐拉法求解m+1步的預(yù)測值：

再用m步微分值f（x（m），y（m））和k+1步預(yù)測的微分值f（x（m+1），y（m+1））的平均值來近似微分：

再用顯式歐拉法校正k+1步的計算值：

綜上所述，改進后的歐拉算法為：

上式就是關(guān)于y（m）向y（m+1）的遞推形式，根據(jù)初始條件按照遞推關(guān)系依次求出y（1），y（2），y（3）…y（n）…，所求出的離散序列為方程的數(shù)值解。

2.4? 三種歐拉方法對比

微分方程數(shù)值解法的本質(zhì)特征是使用數(shù)值積分實現(xiàn)向y的轉(zhuǎn)換。顯式歐拉法和隱式歐拉法都使用的是矩形數(shù)值積分方法，是一階算法，截斷誤差為O（h2）;改進歐拉法使用的是梯形數(shù)值積分方法，是二階算法，截斷誤差為O（h3）。

3? 四階龍格庫塔法

歐拉法在求解方程時，隨著步長的增加，歐拉法一階算法精度下降很快。而龍格庫塔法都是由合適的泰勒方法推導(dǎo)得到，使得誤差為O（h5），其求解精度高，且易于編程。

四階龍格庫塔法的求解算法為式（10）所示。

四階龍格庫塔法的算法是關(guān)于y（m）向y（m+1）的遞推形式，可以根據(jù)初始條件按照遞推關(guān)系依次求出y（1），y（2），y（3）…y（n）…，則離散序列就是方程的數(shù)值解。

4? 數(shù)值實驗

4.1? MATLAB中微分方程數(shù)值中的函數(shù)

在求解微分方程的數(shù)值解時，ode是MATLAB軟件中專有的求解微分方程的函數(shù)。在MATLAB中共有7個函數(shù)，分別是：ode45，ode23，ode113，ode15s，ode23s，ode23t和ode23tb。其中ode45采用了四階-五階Runge-Kutta算法，用四階龍格庫塔法提供候選解，五階龍格庫塔法控制誤差。

4.2? 數(shù)值舉例

例1，求解一階常微分方程：

解：使用MATLAB編程計算得到不同時間步長h的數(shù)值計算的結(jié)果分別如圖1、圖2、表1和表2所示。

由例題數(shù)據(jù)分析可得，顯式歐拉法計算量小，但精度低;而隱式歐拉法的計算遠比顯式歐拉法困難，但其計算穩(wěn)定性好。改進歐拉法的精度高于前兩者。

四階龍格庫塔法即使在步長h很大時，求解的精準(zhǔn)度比MATLAB自帶的ode45函數(shù)和改進歐拉法精度明顯提高。四階龍格庫塔法相對于使用ode45函數(shù)最大的優(yōu)勢在于：可以將求解程序和模型描述文件融合起來，能夠解決各類參數(shù)時變、多模型切換等問題。

5? 結(jié)? 論

通過上述例子可以直觀地看到各種方法在解微分方程數(shù)值解時的優(yōu)缺點，無論是顯式、隱式還是改進歐拉法，精準(zhǔn)度主要取決于區(qū)間的大小和步長h，當(dāng)取得區(qū)間過大時，會出現(xiàn)前一部分區(qū)間的近似值偏差小，后一部分的結(jié)果偏差值較大。

四階龍格庫塔法比顯式歐拉法、隱式歐拉法以及改進歐拉法的精準(zhǔn)度更高，與MATLAB軟件中的ode45函數(shù)相比有著更多的優(yōu)勢。

通過以上對比分析，在求解微分方程的數(shù)值解時，可以供讀者快速便捷的選取合理的算法應(yīng)用于計算，使微分方程在對實際問題的研究應(yīng)用中更加高效快捷。

參考文獻：

[1] 姜兆檸.常微分方程數(shù)值解的求解 [J].科技經(jīng)濟導(dǎo)刊，2019，27（17）：154-155.

[2] 劉相國，楊曉偉，王冬銀.基于MATLAB的《常微分方程》教學(xué)研究 [J].西安文理學(xué)院學(xué)報（自然科學(xué)版），2020，23（2）：124-128.

[3] 王高雄，周之銘.常微分方程：第3版 [M].北京：高等教育出版社，2006.

[4] 梁春葉，王橋明，孫遠通，等.微分方程數(shù)值解之歐拉法在MATLAB下的應(yīng)用 [J].科技風(fēng)，2021（10）：71-72.

[5] 賈新輝.基于微分方程模型信賴域子問題的改進歐拉法 [D].太原：太原科技大學(xué)，2017.

作者簡介：任春年（1999.12—），男，漢族，陜西榆林人，本科在讀，研究方向：數(shù)學(xué)與應(yīng)用數(shù)學(xué)。

3780500338226