龍格庫塔方法是數值求解常微分非線性方程的有利工具,計算精度較高,通過縮短步進距離和增加階數可以進一步控制誤差範圍。工程上較為常用的是四階龍格庫塔演算法(r-k4),在計算收斂的情況下往往可以得到比較好的結果。
這裡簡單介紹一下演算法的具體實現過程,不做詳細的推導。其求解的問題是形如方程:
y ˙=
f(y,
t),其
中t∈[
t0,t
1]初值
y(t0
)=c0
\dot=f(y,t),其中t\in[t_0,t_1] \\ 初值 y(t_0)=c_0
y˙=f(
y,t)
,其中t
∈[t0
,t1
]初值
y(t0
)=c
0通過選取一定的步進長度h,來對區間上函式值進行單步迭代求解,最終得到結果。具體計算公式為:
t n+
1=tn
+hk1
=f(y
n,tn
)k2=
f(yn
+h2k
1,tn
+h2)
k3=f
(yn+
h2k2
,tn+
h2)k
4=f(
yn+h
k3,t
n+h)
yn+1
=yn+
h6(k
1+2k
2+2k
3+k4
)t_=t_n+h\\ k_1=f(y_n,t_n)\\ k_2=f(y_n+\dfrack_1,t_n+\dfrac)\\ k_3=f(y_n+\dfrack_2,t_n+\dfrac)\\ k_4=f(y_n+hk_3,t_n+h)\\ y_=y_n+\frac(k_1+2k_2+2k_3+k_4)
tn+1=
tn+
hk1
=f(y
n,t
n)k
2=f
(yn
+2h
k1,
tn+
2h)
k3=
f(yn
+2h
k2
,tn
+2h
)k4
=f(y
n+h
k3,
tn+
h)yn
+1=
yn+
6h(
k1+
2k2
+2k3
+k4
)通過以上公式,選取合適的步進長度h,反覆迭代,就可求解出y的數值解。
clear
;clc;
close all;
h=1e-5; %步進長度
t=0:h:1; %生成自變數t的向量
%%建立計算結果x,y,z的陣列
n=length(t)
;x=ones(1,n)
;y=ones(1,n)
;z=ones(1,n)
;%%四階龍格庫塔迭代
for i=2:n
t_n=t(i-1)
; x_n=x(i-1)
; y_n=y(i-1)
; z_n=z(i-1)
;
kx1=y_n+3*z_n+sin(5*t_n)
; ky1=x_n+cos(t_n)
; kz1=x_n+z_n-3*cos(3*t_n)*sin(4*t_n)
;
kx2=
(y_n+ky1*h/2)+3*(z_n+kz1*h/2)+sin(5*(t_n+h/2))
; ky2=
(x_n+kx1*h/2)+cos(t_n+h/2)
; kz2=
(x_n+kx1*h/2)+(z_n+kz1*h/2)-3*cos(3*(t_n+h/2))*sin(4*(t_n+h/2))
;
kx3=
(y_n+ky2*h/2)+3*(z_n+kz2*h/2)+sin(5*(t_n+h/2))
; ky3=
(x_n+kx2*h/2)+cos(t_n+h/2)
; kz3=
(x_n+kx2*h/2)+(z_n+kz2*h/2)-3*cos(3*(t_n+h/2))*sin(4*(t_n+h/2))
;
kx4=
(y_n+ky3*h)+3*(z_n+kz3*h)+sin(5*(t_n+h))
; ky4=
(x_n+kx3*h)+cos(t_n+h)
; kz4=
(x_n+kx3*h)+(z_n+kz3*h)-3*cos(3*(t_n+h))*sin(4*(t_n+h))
;
x(i)
=x_n+h/6*(kx1+2*kx2+2*kx3+kx4)
; y(i)
=y_n+h/6*(ky1+2*ky2+2*ky3+ky4)
; z(i)
=z_n+h/6*(kz1+2*kz2+2*kz3+kz4)
; end
%%畫圖
figure();
hold on;
plot(t,x,'r');
plot(t,y,'g');
plot(t,z,'b');
legend(
'x','y','z');
xlabel(
't')
;title(
'xyz函式圖象');
hold off;
四階龍格 庫塔(Runge Kutta)方法
簡潔版 code fprintf 請輸入區間下界 n a input fprintf 請輸入區間上界 n b input fprintf 請輸入初值alpha n alpha input fprintf 請輸入最大迭代次數n n n input x0 a y0 alpha h b a n k x1 ...
四階龍格庫塔法
這裡主要講一下如何用c語言程式設計運用四階龍格庫塔法求解微分方程組。對於所舉例子,只是為了說明龍格庫塔法不僅可以解一階線性微分方程,高階非線性也可通過降階後按照經典四階龍格庫塔法公式逐步求解。只要選取合適的步長h,就能夠平衡速度和精度,達到求解要求。至於例子中的一級倒立擺的物理含義沒有提及到,各種方...
四階龍格庫塔法的基本思想 四階龍格庫塔實驗報告
1 三 四階runge kutta法求解常微分方程 一 龍格庫塔法的思想根據第九章的知識可知道,euler方法的區域性截斷誤差是,而當用euler方法估計出再用梯形公式進行校正,即採用改進euler方法得出數值解的截斷誤差為。由lagrange微分中值定理記,得到這樣只要給出一種計算的演算法,就能得...