四階龍格庫塔法

這裡主要講一下如何用c語言程式設計運用四階龍格庫塔法求解微分方程組。對於所舉例子，只是為了說明龍格庫塔法不僅可以解一階線性微分方程，高階非線性也可通過降階後按照經典四階龍格庫塔法公式逐步求解。只要選取合適的步長h，就能夠平衡速度和精度，達到求解要求。至於例子中的一級倒立擺的物理含義沒有提及到，各種方程具體是如何來的也沒有涉及到推導。關於控制理論方面的知識這裡不作重點。

對微分方程：dy/dt = f(x, y)

有初值條件：y（x(i)）= φ（x(i)）

y(i+1) = y(i)+h*(k1+2*k2+2*k3+k4)/6

k1=f(x(i),y(i))

k2=f(x(i)+h/2,y(i)+h*k1/2)

k3=f(x(i)+h/2,y(i)+h*k2/2)

k4=f(x(i)+h,y(i)+h*k3)

其中，k1， k2， k3， k4表示的是輸出變數的一階倒數，即在一點處的微分，斜率。

以下以乙個直線一級倒立擺為例。

根據牛頓經典力學可以建立以下二階非線性微分方程組模型：

其狀態方程為：

其輸出方程為：

由微分方程可知，是乙個高階微分方程，先對其進行降階：

令不管是c語言程式設計還是matlab程式設計，都只需根據公式和所要求解的具體微分方程來求解各斜率，迴圈4次，最終求解輸出變數。

for(i = 0;i < 4;i++)			
if(i==2)
}

最終計算x值程式語句：

x[0] = x[0]+h*(k[0][0]+2*k[0][1]+2*k[0][3]+k[0][3])/6;
x[1] = x[1]+h*(k[1][0]+2*k[1][1]+2*k[1][3]+k[1][3])/6;
x[2] = x[2]+h*(k[2][0]+2*k[2][1]+2*k[2][3]+k[2][3])/6;
x[3] = x[3]+h*(k[3][0]+2*k[3][1]+2*k[3][3]+k[3][3])/6;

注意選取合適的步長以及計算區間。

四階龍格庫塔法

四階龍格庫塔法的基本思想四階龍格庫塔實驗報告

python實現四階龍格庫塔法

四階龍格庫塔（Runge Kutta）方法

四階龍格庫塔法

四階龍格庫塔法的基本思想 四階龍格庫塔實驗報告

python實現四階龍格庫塔法

四階龍格 庫塔（Runge Kutta）方法

相關推薦

四階龍格庫塔法的基本思想四階龍格庫塔實驗報告

四階龍格庫塔（Runge Kutta）方法