四階龍格庫塔法

這里主要講一下如何用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++)

{

//直接寫出所要求解的微分方程

dX[0] = x2;

dX[1] = u;

dX[2] = x4;

dX[3] = (3*u*cos(x3))+(29.4*sin(x3))-(0.165*(x4));

//將倒數即斜率賦值給k

k[0][i] = DX[0];

k[1][i] = DX[1];

k[2][i] = DX[2];

k[3][i] = DX[3];

//計算出x的值，供求下一組斜率k時使用

if(i==0||i==1)

{

x1 = X[0]+h*k[0][i]/2;

x2 = X[1]+h*k[1][i]/2;

x3 = X[2]+h*k[2][i]/2;

x4 = X[3]+h*k[3][i]/2;

}

if(i==2)

{

x1=X[0]+h*k[0][i];

x2=X[1]+h*k[1][i];

x3=X[2]+h*k[2][i];

x4=X[3]+h*k[3][i];

}

}

最終計算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; 注意選取合適的步長以及計算區間。

總結

以上是生活随笔為你收集整理的四阶龙格库塔法程序c语言,四阶龙格库塔法的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。