Runga-Kutta Method for system of first order differential equations

68 views (last 30 days)
I am trying to work out the runga kutta method for a system of first order ODE's. I have the scheme already set up and it runs without any errors, giving an estimation. The ODE's I am trying to solve:
This scheme blows up, which is not what I was expecting. Is there an error in the scheme? Or have I missunderstood the fact that this scheme should not blow up?
Below you can find the code I have written:
%% Initialization
tspan = 100;
h = 0.125;
fact = 1/h;
x = 0:h:tspan;
y = zeros(length(x),2);
y(1,1) = 0.2; % Initial condition y1(0) = 0.2
y(1,2) = 0; % Initial coniditon y2(0) = 0
dy1dt = @(t,y) y(:,2); % change the function as you desire
dy2dt = @(t,y) -y(:,1); % change the function as you desire
%% Runge Kutta
if Method == 'RKM' | Method == 'All'
for i3 = 1:(length(x)-1) % calculation loop
k_1 = dy1dt(x(i3) , y(i3,:));
k_2 = dy1dt(x(i3) + 0.5*h , y(i3,:) - 0.5.*h.*k_1);
k_3 = dy1dt((x(i3) + 0.5*h) , (y(i3,:) - 0.5.*h.*k_2));
k_4 = dy1dt((x(i3) + h) , (y(i3,:) - k_3.*h));
y(i3+1,1) = y(i3,1) + (1/6)*(k_1 + 2*k_2 + 2*k_3 + k_4).*h; % main equation
l_1 = dy2dt(x(i3) , y(i3,:));
l_2 = dy2dt(x(i3) + 0.5*h , y(i3,:) - 0.5*h*l_1);
l_3 = dy2dt((x(i3) + 0.5*h) , (y(i3,:) - 0.5*h*l_2));
l_4 = dy2dt((x(i3) + h) , (y(i3,:) - l_3*h));
y(i3+1,2) = y(i3,2) + (1/6)*(l_1 + 2*l_2* + 2*l_3 + l_4)*h; % main equation
end
figure(3),clf(3),hold on
plot(x,y(:,1),'r')
plot(x,y(:,2),'b')
plot(T,Y(:,1),'r--') % Solution according to ODE45
plot(T,Y(:,2),'b--') % Solution according to ODE45
axis([0 100 -1 1])
legend('Runga Kutta 1','Runga Kutta 2','Exact 1','Exact 2')
title('Runga Kutta')
grid on;
end
The result of this code is below.

Accepted Answer

Davide Masiello
Davide Masiello on 5 May 2022
Edited: Davide Masiello on 5 May 2022
Two things:
1 - By defining two different functions for y1 and y2 you decouple the system, which is the source of the main problem.
2 - After fixing point 1, I had to use a very small time step size to obtain an acceptable solution, indicating that the implementation of an adaptive step size might be beneficial to the efficiency of the code.
%% Initialization
tspan = 100;
h = 0.0005;
fact = 1/h;
x = 0:h:tspan;
y = zeros(2,length(x));
y(:,1) = [0.2;0]; % Initial conditions
%% Runge Kutta
for i = 1:(length(x)-1) % calculation loop
k_1 = odeSystem(x(i),y(:,i));
k_2 = odeSystem(x(i)+0.5*h,y(:,i)-0.5.*h.*k_1);
k_3 = odeSystem((x(i)+0.5*h),(y(:,i)-0.5*h*k_2));
k_4 = odeSystem((x(i)+h),(y(:,i)-k_3*h));
y(:,i+1) = y(:,i)+(1/6)*(k_1 + 2*k_2 + 2*k_3 + k_4)*h; % main equation
end
figure
plot(x,y(1,:),'r',x,y(2,:),'b')
axis([0 100 -1 1])
title('Runge-Kutta')
function dydt = odeSystem(t,y)
dydt(1,1) = y(2);
dydt(2,1) = -y(1);
end

More Answers (1)

Chunru
Chunru on 5 May 2022
Looks like you need to choose a much smaller step size for the problem.
%% Initialization
tspan = 100;
h = 0.125/50;
fact = 1/h;
x = 0:h:tspan;
y = zeros(length(x),2);
y(1,1) = 0.2; % Initial condition y1(0) = 0.2
y(1,2) = 0; % Initial coniditon y2(0) = 0
dy1dt = @(t,y) y(:,2); % change the function as you desire
dy2dt = @(t,y) -y(:,1); % change the function as you desire
%% Runge Kutta
Method = "RKM"
Method = "RKM"
if Method == 'RKM' | Method == 'All'
for i3 = 1:(length(x)-1) % calculation loop
k_1 = dy1dt(x(i3) , y(i3,:));
k_2 = dy1dt(x(i3) + 0.5*h , y(i3,:) - 0.5.*h.*k_1);
k_3 = dy1dt((x(i3) + 0.5*h) , (y(i3,:) - 0.5.*h.*k_2));
k_4 = dy1dt((x(i3) + h) , (y(i3,:) - k_3.*h));
y(i3+1,1) = y(i3,1) + (1/6)*(k_1 + 2*k_2 + 2*k_3 + k_4).*h; % main equation
l_1 = dy2dt(x(i3) , y(i3,:));
l_2 = dy2dt(x(i3) + 0.5*h , y(i3,:) - 0.5*h*l_1);
l_3 = dy2dt((x(i3) + 0.5*h) , (y(i3,:) - 0.5*h*l_2));
l_4 = dy2dt((x(i3) + h) , (y(i3,:) - l_3*h));
y(i3+1,2) = y(i3,2) + (1/6)*(l_1 + 2*l_2* + 2*l_3 + l_4)*h; % main equation
end
figure(3),clf(3),hold on
plot(x,y(:,1),'r')
plot(x,y(:,2),'b')
% plot(T,Y(:,1),'r--') % Solution according to ODE45
% plot(T,Y(:,2),'b--') % Solution according to ODE45
axis([0 100 -1 1])
legend('Runga Kutta 1','Runga Kutta 2','Exact 1','Exact 2')
title('Runga Kutta')
grid on;
end
Warning: Ignoring extra legend entries.

Products


Release

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!