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

梯形積分法

f(x)在封閉區間[a,b]求積分，若f(x)為複雜的函數， 對於求解積分，可用

數值分析的梯形積分法來求解。

梯形積分法，得到的結果，對照原函數定積分的值，可發現誤差值很小。所以被廣泛運用。

以下用matlab 語法，撰寫梯形積分的algorithm：

%Function F(x)=∫sin(x)dx，積分範圍0~π

a=0;

b=pi;

n=1000;

A=0;

for i=a:(b-a)/n:b

　　if i==a

　　　f1=sin(i);

　　else

　　　f2=sin(i);

　　　A=A+0.5*(f1+f2)*(b-a)/n;

　　　f1=f2;

　　end

end

disp('梯形積分法面積=');

A

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

經執行後，A＝2.0000

與積分範圍[0,pi]，∫sin(x)dx=2的結果一致。

Runge Kutta Method order 2

微分方程式y'=y； y(0)=1

Ans:

使用Runge Kutta Method order 2:

Runge-Kutta Method通式：

y'＝f(x,y)

x_i=x₀+ih,i=0,1,2…

y_k+1＝y_k＋(1/2)．(k₁+k₂)

其中

k₁=h．f(x_k,y_k)

k₂=h．f(x_k＋h,y_k＋k1)

假設h=0.5(漸進步距)

所以：

y'＝f(x,y)＝y

y(0)＝1，換句話說y₀＝1

x₁=x₀+i×h=0+1×0.5=0.5

k₁＝h×f(x₀,y₀)＝0.5×f(0,1)=0.5×1=0.5

k₂＝h×f(x₀＋h,y₀＋k₁)＝0.5×f(0＋0.5,1＋0.5)=0.5×f(0.5,1.5)=0.5×1.5=0.75

其中f(x,y)=y，帶入x=0,y=1進f(x,y)中，得f(0,1)=1

--> y₁＝y₀＋(1/2)×(k₁+k₂)=1+(1/2)×(0.5+0.75)=1.625

x₂=x₀+i×h=0+2×0.5=1.0

k₁＝h×f(x₁,y₁)＝0.5×f(0.5,1.625)=0.5×1.625=0.8125

k₂＝h×f(x₁＋h,y₁＋k₁)

＝0.5×f(0.5＋0.5,1.625＋0.8125)

＝0.5×f(1,2.4375)=0.5×2.4375=1.2188

-->y₂＝y₁＋(1/2)×(k₁+k₂)=1.625+(1/2)×(0.8125+1.2188)=2.6406

對正式解而言，前述y’=y，y(0)=1的通解是y=e^x：

1. y(0)=1

2. y(0.5)=1.6487，經由Runge-Kutta方法計算的y1=1.625，

誤差率＝(1.6487-1.625)/1.625=1.46%

3.y(1)=2.71828，經由Runge-Kutta方法計算的y2=2.6406，

誤差率＝(2.71828-2.6406)/2.6406=2.94%