Runge-kutta 4

Please I need your help, I was told to solve this with euler forward method and Runge kutta method. Sadly i could only use Euler forward method, please how can this be done using Runge kutta method, sorry i am writing here, i have a deadline for Friday midnight. Thanks

#Edit: formatted your code

b=0.0020 

b = 0.0020

gama=0.02; 
D=30; 
dt=1/60; 
nt=floor((D*24)/dt) 

nt = 43200

t=linspace(0,nt*dt,nt+1); 
S=zeros(nt+1,1); 
I=zeros(nt+1,1); 
R=zeros(nt+1,1); % this is the initial conditiion 
S(1)=50; %MATLAB don't recognise 0 
I(1)=1; 
R(1)=1; 
for i=1:nt 
    S(i+1)=S(i)-dt*b*S(i)*I(i); 
    I(i+1)=I(i)+dt*b*S(i)*I(i)-dt*gama*I(i); 
    R(i+1)=R(i)+dt*gama*I(i); 
end 
figure(1) 
plot(t,S,t,I,t,R,'LineWidth',2); 
legend('S','I','R'); 
xlabel('hours'); 
ylabel('population');

[Sam Chak]

Can you show the Runge–Kutta 4 iterative algorithm in math form?

[Stanley Umeh]

Sorry I do not know what you mean, I am trying to do it with the rungs kutta, and sadly It's no working, so I watching video, can't this euler method help you with the runge kutta?

[Jan]

@Stanley Umeh: Please post your current code and explain "not working" with any details. Then the readers can suggest a modification of the code.

You can find dozens of Runge-Kutta implementations in Matlab. Simply search in the FileExchange or in the net.

The forum will help you, but it will not solve your homework.

[Stanley Umeh]

Writing k1 and k2 is my problem, how can I arrange k1 and k2 using runger kutta method.

b=0.0015; gama=0.01; D=30; %Day dt=1/60; %time increment (Minutes) nt=floor((D*24)/dt); t=linspace(0,nt*dt,nt+1); S=zeros(nt+1,1); I=zeros(nt+1,1); R=zeros(nt+1,1); % initial conditions S(1)=50; I(1)=1; R(1)=1; for i=1:nt S(i+1)=S(i)-dt*b*S(i)*I(i); K1= K2=

[Jan]

The code will be easier, if you join S, I, R to a vector.

Remember, that there are a lot of different Runge-Kutta methods. You are asking for a method of order 4. Then "k1 and k2" are not enough.

Did you follow my advice and searched in the FileExchange already? MathWorks has published fixed step solvers also: https://www.mathworks.com/matlabcentral/answers/98293-is-there-a-fixed-step-ordinary-differential-equation-ode-solver-in-matlab-8-0-r2012b#answer_107643

[Stanley Umeh]

With the help of k1 and k2 i can do k3 and 4. I will check it our, my code for euler is complete just need help with RK 4th order. I will check it

Answer 1

[Sam Chak]

Hi @Stanley Umeh

I requested you to show the mathematics of RK4 so that you can learn about how the algorithm works as well as saving time for me to do a google search and type out the Rk4 formulas. The RK4 code may contain errors. Can you check if the formulas in the code are correct?

% Parameters
b    = 0.0020;
gama = 0.02; 

% SIR model (Ordinary Differential Equations)
odefcn = @(t,x) [- b*x(1)*x(2);
                 b*x(1)*x(2) - gama*x(2);
                 gama*x(2)];

% Simulation time span
tStart = 0;
h      = 0.2;   % step size
tEnd   = 800;
t      = tStart:h:tEnd;	        

% Initial values
y0     = [50 1 1];

% Calling RK4 Solver is similar to ode45 solver
yRK4 = rk4(odefcn, t, y0);

% Plotting the solutions
plot(t, yRK4, 'linewidth', 1.5)
grid on
xlabel('Hours')
ylabel('Population')
title('Time responses of SIR model')
legend('S(t)', 'I(t)', 'R(t)', 'location', 'best')

% Code for Runge-Kutta 4th-order Solver
function y = rk4(f, x, y0)
    y(:, 1) = y0;               % initial condition
    h       = x(2) - x(1);      % step size
    n       = length(x);        % number of steps
        for j  = 1 : n-1
            k1 = f(x(j), y(:, j)) ;
            k2 = f(x(j) + h/2, y(:, j) + h/2*k1) ;
            k3 = f(x(j) + h/2, y(:, j) + h/2*k2) ;
            k4 = f(x(j) + h, y(:, j) + h*k3) ;
            y(:, j+1) = y(:, j) + h/6*(k1 + 2*k2 + 2*k3 + k4) ;
        end
end