Runge Kutta Method order 4-Example 1

Example 1: solve dy/dx=f(x,y)=3exp(-x)-0.4y，y(0)=5，find y(3)=？

Ans:

使用Runge Kutta Order 4求近似解，已知Runge Kutta Method order 4的一般式：

dy/dx=f(x,y)

→y(k+1)=y(k)＋(1/6)×h×(k1+2k2+2k3+k4)

其中:

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

k2=f(x(k)+(1/2)h,y(k)+(1/2)k1×h)

k3=f(x(k)+(1/2)h,y(k)+(1/2)k2×h)

k4=f(x(k)+h,y(k)+k3×h)

h: 遞增步距

程式範例：

%dy/dx=f(x,y)=3exp(-x)-0.4y, y(0)=5

%Using Runge-Kutta order 4: y(k+1)=y(k)+1/6(k1+2k2+2k3+k4)h

clear all;

clc

h=0.1; %Step

L=0;

R=3;

X=0;Y=5;

counter=1;

for v=L+h:h:R

　k1=Func(X,Y);

　k2=Func(X+0.5*h,Y+0.5*k1*h);

　k3=Func(X+0.5*h,Y+0.5*k2*h);

　k4=Func(X+h,Y+k3*h);

　Y=Y+h*(k1+2*k2+2*k3+k4)/6;

　Y_sol=10*exp(-0.4*(X+h))-5*exp(-(X+h));

　err(counter)=Y_sol-Y;

　Yk(counter)=Y;

　Yk_sol(counter)=Y_sol;

　counter=counter+1;

　X=X+h;

　plot(Yk,'r-')

　hold on

　plot(Yk_sol,'b*')

end

hold off

副程式func.m

function k=Func(x,y)

k=3*exp(-x)-0.4*y;

end

----------------------------------

執行結果：

h=0.1（分析步距）

原微分方程式的解y(x)=10e^-0.4x-5e^-x

上圖中，紅色線為期望值，藍色點是近似值。

利用RK4分析微分方程式的動態行為，其近似值逼近期望值。

總誤差=Σ(每一點的誤差)=1.3587×10^-6