四階龍格庫塔方法求解一次常微分方程組

2021-10-16 22:40:41 字數 2872 閱讀 7334

龍格庫塔方法是數值求解常微分非線性方程的有利工具,計算精度較高,通過縮短步進距離和增加階數可以進一步控制誤差範圍。工程上較為常用的是四階龍格庫塔演算法(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微分中值定理記,得到這樣只要給出一種計算的演算法,就能得...