科學計算—理論、方法
及其基於matlab的程式實現與分析
微分方程（組）數值解法
§1 常微分方程初值問題的數值解法
微分方程（組）是科學研究和工程應用中最常用的數學模型之一。如揭示質點運動規律的newton第二定律：
和刻畫迴路電流或電壓變化規律的基爾霍夫迴路定律等，但是，只有一些簡單的和特殊的常微分方程及常微分方程組，可以求得用公式給出的所謂「解析解」或「公式解」，如一階線性微分方程的初值問題：
的解為：
但是，絕大多數在實際中遇到的常微分方程和常微分方程組得不到「解析解」，因此，基於如下的事實：
１、絕大多數的常微分方程和常微分方程組得不到（有限形式的）解析解；
２、實際應用中往往只需要知道常微分方程（組）的解在（人們所關心的）某些點處的函式值（可以是滿足一定精度要求的近似值）；
如果只需要常微分方程（組）的解在某些點處的函式值，則沒有必要非得通過求得公式解，然後再計算出函式值不可，事實上，我們可以採用下面將介紹的常微分方程（組）的初值問題的數值解法，就可以達到這一目的。
一般的一階常微分方程（組）的初值問題是指如下的一階常微分方程（組）的定解問題：
其中常微分方程（組）的初值問題通常是對一動態過程（動態系統、動力系統）演化規律的描述，求解常微分方程（組）的初值問題就是要了解和掌握動態過程演化規律。
§１.1 常微分方程（組）的cauch問題數值解法概論
假設要求在點（時刻），處初值問題（７）的解的（近似）值,如果已求得時刻的值或它的近似值(如時刻的值)，那麼將式（７）的兩端在區間上積分
(10)
可得11)
或 (12)
顯然,為了利用式（11）或（12）求得的精確值（近似值），必須計算右端的積分，這是問題的關鍵也是難點所在，如前所述，一般得不到精確的公式解，因此需要採用數值積分的方法求其近似解,
可以說,不同的式值積分方法將給出不同的cauch問題的數值解法。
§1.2 最簡單的數值解法——euler 方法
假設要求在點（時刻），，處初值問題（７）的解的近似值。首先對式（７）的兩端積分，得
（13）
對於式（13）的右邊，如果用積分下限處的函式值代替被積函式作積分（從幾何上的角度看，是用矩形面積代替曲邊梯形面
積），則有
（14）
進而得到下式給出的遞推演算法—euler 方法
（15）
例１ 用euler 方法解如下初值問題，取，
解：由（15）得
結果如下：
open euler_method.m
如果取，其結果如下圖所示：
euler_method
§1.3 改進的euler 方法
對於（15）的右邊，如果被積函式用積分限和處的函式值的算術平均值代替（幾何上，是用梯形面積代替曲邊梯形面積），則有
（16）
進而得到下式給出的遞推演算法：
（17）
通常演算法（17）比euler 方法（15）的精度高，但是，按演算法（17）求時要解（非線性）方程（組），這是演算法（17）不如euler 方法的方面，為了
1） 盡可能地保持演算法（17）精度高的優點；
2） 盡可能地利用euler 方法計算簡單的長處；
人們採取了如下的稱之為改進的euler 方法的折衷方案：
**18）
修正 （19）
例２ euler 方法與改進的euler 方法的比較
下圖是當時比較的結果：
open improved_euler_method.m
§1.4 euler 方法和改進的euler 方法的誤差分析
由taylor 公式
（19）
說明euler 方法的截斷誤差是，類似地，由
20）21）
以及22)
讓式（20）的兩端減式（21）的兩端，可得
（23）
從上述推導euler 方法、改進的euler 方法的過程以及例１、例２容易看出，改進的euler 方法euler 方法的精度高，其原因在於：
1 在推導euler 方法時，我們是用待求解函式在一點處的變化
率代替在區間上的平均變化率：
24）2 而在推導改進的euler 方法時，我們是用待求解函式在兩點處變化率的平均值代替在區間上的平均變化率；
顯然，通常比更接近於在區間上的平均變化率。由此啟發人們：適當地選取區間上函式若干點處的變化率，用它們加權平均值代替在區間上的平均變化率，近似解的精度應更高。
下面將要介紹的runge—kutta法就是基於上述想法得到的。
§2 runge—kutta法
runge—kutta法是按選取區間上函式變化率的個數的多少和截斷誤差的階數來區分的一系列方法，如
１ 二階的runge—kutta法（改進的euler 方法）
25）２ 三階的runge—kutta法
26）３ 四階的runge—kutta法
1) 古典形式
27） 2) gill 公式(具有減小捨入誤差的優點)
28）4 runge—kutta法的一般形式
29)其中稱為增量函式(increment function)以及
（30）
需要特別指出的是：在確定（29）、（30）中的引數，和時，應該使（29）的右端和適當階的（如階）taylor展式一致，這樣，至少對於底階的r-k法來說，參與加權的斜率個數與方法的階數是一致的。
例如：一階的r-k法即euler法，區域性的截斷誤差（truncation error）是，所以，當微分方程的階是一次函式時是精確的；二階r-k法即改進的euler法，區域性的截斷誤差是，所以，當微分方程的階是二次函式時是精確的；類似地，四階r-k法區域性的截斷誤差是。
但是，一般五階以及五階以上的r-k法的區域性截斷誤差不再具有上述的特點，因此用「價效比」來衡量，四階r-k得到了廣泛的應用。
下面是butcher's fifth-order r-k 方法（1964）：
28） （29）
例三：疾病的傳播與防疫模型
open lunge_kutta_method.m
其中表示一種傳染病傳染者的人數，表示易感染者人數。
31)該系統的jacobi矩陣為：
32)該系統的平衡點為：
33)根據問題的實際意義，有意義的平衡點為：
34)§3 解剛性問題的數值方法
例四：單個微分方程的例子
考慮如下的一階微分方程描述的系統
35)容易求的該微分方程初值問題的通解：
36)f=dsolve('dy = -1000*y+3000-2000*exp(-t)', 'y(0) = c')
並且有:
37)當初值時，方程的解中包含「快變」和「慢變」兩部分.
應用不同數值方法求解的結果：
１ 解非剛性問題的方法：ode45sol_stiffode01.m
２ 解剛性問題的方法： ode15ssol_stiffode01.m
３ 解非剛性問題的方法：euler method: sol_stiffode01_euler.m
系統（35）所對應的自治系統為
38)利用euler 法求解上述初值問題：
39)為保證時間序列的收斂性，步長應滿足不等式：.
例五：微分方程組的例子
考慮如下的一階微分方程組描述的系統
40) (41)
(42)
利用euler 法求解上述初值問題的迭代公式為
43)由於係數矩陣的特徵值為
(44)
所以euler 法不適用於該初值問題的求解。
利用四階r-k 法求解上述初值問題的迭代公式為
45)其中係數矩陣為
(46)
其特徵值為：
(47)
１ 解非剛性問題的方法：ode45sol_stiffode02.m
２ 解剛性問題的方法： ode15ssol_stiffode02.m
３ 解非剛性問題的方法：euler method: sol_stiffode02_euler.m
４ 解非剛性問題的方法：euler method: sol_stiffode02_rk4.m
f=dsolve('dx=-5*x+3*y','dy=100*x-301*y','x(0)=1','y(0)=2')
transient
dominate
motivation
orientation
scope
discretization
predictor – corrector approach
increment function
recurrence
adaptive steps - size control
abrupt change
例2：二體問題
open runge_kutta_method01
§4 線性多步法與常微分方程數值解法的分類
一般的常微分方程初值問題的數值解法都是以遞推的形式給出的，即是遞推演算法，根據遞推演算法，在計算時，已經得到了前面各時刻的近似值，，前面介紹的各種數值解法都有乙個共同點：在計算時，只用到了前乙個時刻（當前時刻）的「資訊」**未來某時刻系統的「狀態」，這樣的數值解法稱為單步法，對於乙個動態過程在時刻的狀態而言，不僅前乙個時刻的資訊對它有影響，前若干個時刻的資訊通常對它也有影響，顯然，單步法的缺點是沒有充分利用已得到的資訊。
計算方法 常微分方程的差分方法實驗
實驗三常微分方程的差分方法實驗 一.實驗目的 1 深入理解常微分方程的差分方法的原理，學會用差分方法解決某些實際的常微分方程問題，比較這些方法解題的不同之處。2 熟悉matlab程式設計環境，利用matlab實現具體的常微分方程。二.實驗要求 用matlab軟體實現尤拉方法 改進的尤拉方法 龍格 庫...
數值分析報告常微分方程的數值解法
常微分方程的數值解法 一 實驗目的 目的與要求 通過實驗，熟悉常微分方程的數值解法的基本原理。掌握向前尤拉法 向後尤拉法 梯形法 改進尤拉法及三階 四階龍格 庫塔法等基本演算法。二 實驗內容 在下列方法中 向前尤拉法 向後尤拉法 梯形法 改進尤拉法及三階 四階龍格 庫塔法 選擇不同的三種演算法求下面...
常微分方程數值解實驗報告
學院 數學與資訊科學 專業 資訊與計算科學 姓名 鄭思義 學號 201216524 課程 常微分方程數值解 實驗一 常微分方程的數值解法 1 分別用euler法 改進的euler法 預報校正格式 和s k法求解初值問題。h 0.1 並與真解作比較。1.1實驗 尤拉法 function x,y nae...