（建議閱讀原文）預備知識中點法解常微分方程（組）龍格庫塔法是一類數值解微分方程的演算法， 其中較常見的是四階龍格庫塔法． 這裡不進行推導， 僅僅給出公式如下（
其中
由以上兩式， 不難把該演算法拓展到方程組的情況． 對於
元微分方程組
我們可以把該式記為向量函式的形式
現在我們僅需要把式 1 和 式 2 中的所有
和 都變為
維列向量
和 即可將微分方程拓展為微分方程組．
例程
到目前為止， 我們每求乙個微分方程的數值解都要重新寫一次程式， 對於一些較為複雜的演算法這樣做效率較低． 我們這裡不妨把四階龍格庫塔法寫到乙個單獨的函式檔案 oderk4.m 中， 當我們要解某個特定的方程時， 只需把式 2 中的
作為自變數輸入即可解出
我們先來看第 1 行的函式宣告， 輸入變數中，f是式 4 中
的函式控制代碼，
tspan是乙個2×1的列向量，tspan(1)是初始時間，tspan(2)是終止時間，y0是乙個列向量，y0(ii)是第ii個因變數的初始值，nt是
的個數，
tspan定義的時間區間被等分為nt - 1個小區間． 因變數中，y的行數是因變數的個數， 列數是nt，t是乙個行向量， 由第 6 行定義，y(ii, jj)就是第ii個變數在t(jj)時刻的值． 第 5 行把初值y0賦給y的第 1 列， 第 8-14 行的迴圈根據式 1 和式 2 的向量形式由y的第ii列（
）求第
ii+1列（
）． 我們先來用這個函式來計算 「天體運動的簡單數值計算」 中的問題． 我們令因變數
的四個分量依次為一階方程組（式 4 ）
中的 ． 程式**如下
執行結果如圖 1 所示．
圖 1：執行結果
我們先來看函式fun（20 行）， 這個函式就相當於式 5 ． 第乙個輸入變數y是乙個列向量， 是
的值， 第二個輸入變數是
， 但由於式 5 中沒有出現
， 我們用波浪線代替． 第三個輸入變數是引數
， 即萬有引力常數和中心天體質量之積． 輸出變數
y1是乙個列向量， 是
的值．
再來看主函式keplerrk4（第 1 行）， 引數設定中除了步數從 4000 變為了 100， 其他都和 「天體運動的簡單數值計算」 中的程式一樣， 然而這裡執行結果卻精確得多（曲線幾乎閉合）， 可見這種演算法的優越性．
主函式第 10 行中將fun(y, t, gm)變為函式控制代碼f(y, t)， 這樣gm就可以在 「引數設定」 中設定， 而不用在fun函式內部設定． 第 11 行呼叫了上文中的oderk4函式解方程組， 由於我們在畫圖時不需要用t， 所以把第二個輸出變數改為波浪線．
二階偏微分方程組 龍格庫塔法 常微分法數值計算方法
其中向量 為狀態變數 構成的向量，即 常稱為系統的狀態向量，n稱為系統的階次，而 為任意函式數，t為時間變數，這樣就可以採用數值方法求解常微分方程組了。另外任意高階微分方程都可以通過變數替換變成一階微分方程組,這裡不再贅述.尤拉法在 時刻系統狀態向量的值為 若選擇計算步長h，則可以寫出在 時刻系統狀...
四階龍格庫塔方法求解一次常微分方程組
龍格庫塔方法是數值求解常微分非線性方程的有利工具，計算精度較高，通過縮短步進距離和增加階數可以進一步控制誤差範圍。工程上較為常用的是四階龍格庫塔演算法 r k4 在計算收斂的情況下往往可以得到比較好的結果。這裡簡單介紹一下演算法的具體實現過程，不做詳細的推導。其求解的問題是形如方程 y f y,t ...
四階龍格庫塔法
這裡主要講一下如何用c語言程式設計運用四階龍格庫塔法求解微分方程組。對於所舉例子，只是為了說明龍格庫塔法不僅可以解一階線性微分方程，高階非線性也可通過降階後按照經典四階龍格庫塔法公式逐步求解。只要選取合適的步長h，就能夠平衡速度和精度，達到求解要求。至於例子中的一級倒立擺的物理含義沒有提及到，各種方...