# Runge-Kutta 4th order method

3,100 views (last 30 days)

Show older comments

Mariam Gasra
on 5 May 2019

% It calculates ODE using Runge-Kutta 4th order method

% Author Ido Schwartz

clc; % Clears the screen

clear;

h=5; % step size

x = 0:h:100; % Calculates upto y(3)

Y = zeros(1,length(x));

y(1) = [-0.5;0.3;0.2];

% initial condition

F_xy = @(t,r) 3.*exp(-t)-0.4*r; % change the function as you desire

for i=1:(length(x)-1) % calculation loop

k_1 = F_xy(x(i),y(i));

k_2 = F_xy(x(i)+0.5*h,y(i)+0.5*h*k_1);

k_3 = F_xy((x(i)+0.5*h),(y(i)+0.5*h*k_2));

k_4 = F_xy((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

display(Y(i+1));

if i run the programme i get answer =0;

how can i solve this problem if i have three initial condition -0.5 ,0.3 and 0.2

while x=0:5:100

and how i can plot the answer with respect to x?

### Accepted Answer

David Wilson
on 6 May 2019

Edited: MathWorks Support Team
on 18 Apr 2023

First up, you will need a much smaller step size to get an accurate solution using this explicit RK4 (with no error control). I suggest h = 0.05. Validate using say ode45 (which does have error control).

Then you will need to run your ode above three separate times, once starting from y(1) = -0.5, again with y(1) = 0.3, etc.

Then finally plot the result with plot(x,y,'o-').

h=0.05; % step size

x = 0:h:100; % Calculates upto y(3)

y = zeros(1,length(x));

%y(1) = [-0.5;0.3;0.2];

y(1) = -0.5; % redo with other choices here.

% initial condition

F_xy = @(t,r) 3.*exp(-t)-0.4*r; % change the function as you desire

for i=1:(length(x)-1) % calculation loop

k_1 = F_xy(x(i),y(i));

k_2 = F_xy(x(i)+0.5*h,y(i)+0.5*h*k_1);

k_3 = F_xy((x(i)+0.5*h),(y(i)+0.5*h*k_2));

k_4 = F_xy((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

% validate using a decent ODE integrator

tspan = [0,100]; y0 = y(1);

[tx, yx] = ode45(F_xy, tspan, y0);

plot(x,y,'o-', tx, yx, '--')

##### 4 Comments

Walter Roberson
on 14 Nov 2022

### More Answers (6)

Sandip Das
on 28 Jul 2021

%Published in 25 July 2021

%Sandip Das

clc;

clear all;

dydt=input('Enter the function : \n');

t0=input('Enter the value of t0 : \n');

y0=input('Enter the value of y0 : \n');

tn=input('Enter the value of t for which you want to find the value of y : \n');

h=input('Enter the step length : \n');

i=0;

while i<tn

k_1 = dydt(t0,y0);

k_2 = dydt(t0+0.5*h,y0+0.5*h*k_1);

k_3 = dydt((t0+0.5*h),(y0+0.5*h*k_2));

k_4 = dydt(((t0)+h),(y0+k_3*h));

nexty = y0 + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;

y0=nexty

t0=t0+h

i=i+h;

end

fprintf('The value of y at t=%f is %f',t0,y0);

##### 0 Comments

mahmoud mohamed abd el kader
on 29 Oct 2020

function [x,y] = rk4th(dydx,xo,xf,yo,h)

x = xo:h:xf ;

y = zeros(1,length(x));

y(1)= yo ;

for i = 1:(length(x)-1)

k_1 = dydx(x(i),y(i));

k_2 = dydx(x(i)+0.5*h,y(i)+0.5*h*k_1);

k_3 = dydx((x(i)+0.5*h),(y(i)+0.5*h*k_2));

k_4 = dydx((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;

end

dydx = @(x,y) 3.*exp(-x)-0.4*y;

%[x,y] = rk4th(dydx,0,100,-0.5,0.5);

%plot(x,y,'o-');

end

##### 3 Comments

soham roy
on 8 Dec 2022

Walter Roberson
on 8 Dec 2022

y = zeros(1,length(x));

would change to

y = zeros(length(x), length(y0));

and below that, each y(INDEX) would be replaced with y(INDEX,:)

Mj
on 7 Nov 2020

Hello everyone!

I have to solve this second order differential equation by using the Runge-Kutta method in matlab:

can anyone help me please? and how can i plot the figure?(a against e)

d2a/de2=(((((2+c2)*(Fu^2))/(1+c2))+1)*(a^c2)-((2+c2/1+c2)*(Fu^2/a))-a^(2+(2*c2)))/(((2+c2)*Fu^2)/(1+c2)*(3+c2));

Fu=1

c2=0 , 0.5 , 1 (there are 3 values for c2)

initial conditions are: a=0.8 , d_a=

David Wilson
on 6 May 2021

Wow, you haven't given us too much to go on, so that makes a real challenge.

First up, your 2nd order ODE is needlessly complex given that Fu=1, and c2 =0 say. (I'm not sure what the other valuesare for, Are you solving this 3 seprate times? (Be good to know if that is the case.)

If you have the symbolic toolbox, it makes it easy to simplify your problem to something doable. First up, I'm going to try and solve it analytically.

syms Fu c2 real

syms a(t)

f2 = (((((2+c2)*(Fu^2))/(1+c2))+1)*(a^c2)-((2+c2/1+c2)*(Fu^2/a))-a^(2+(2*c2)))/(((2+c2)*Fu^2)/(1+c2)*(3+c2));

f2_a = subs(f2,Fu,1)

f2_b = subs(f2_a,c2,0) % subs c2 for 0

Da = diff(a);

D2a = diff(a,2);

% Now attempt to solve analytically

dsolve(D2a == f2_b, a(0) == 0.8, Da(0) == 1)

Well that didn't work, but no real suprise there.

Let's try a numerical method:

syms Fu c2 real

syms a real

f2 = (((((2+c2)*(Fu^2))/(1+c2))+1)*(a^c2)-((2+c2/1+c2)*(Fu^2/a))-a^(2+(2*c2)))/(((2+c2)*Fu^2)/(1+c2)*(3+c2));

f2_a = subs(f2,Fu,1); f2_b = subs(f2_a,c2,0); pretty(f2_b)

We need to encode this as a system of 2 ODES. (Convert to Cauchy form)

aprime = @(t,a) [a(2); ...

0.5 - a(1).^2/6 - 1./(a(1)*3)]

Now we are ready to solve the ODE. I'll use ode45, and guess a t-span, and guess one of the initial conditions since you forgot to help us out there.

aprime = @(t,a) [a(2); ...

0.5 - a(1).^2/6 - 1./(a(1)*3)]

a0 = [0.8; 0]

[t,a] = ode45(aprime, [0,4], a0)

plot(t,a)

##### 0 Comments

Amr Mohamed
on 9 May 2021

how can we write the code for this problem :

##### 2 Comments

Moneeb Ur Rehman
on 27 May 2021

monsef
on 17 Jul 2023

y=x^2-2yx

h=0.2

y0=0

x0=1

wriet program im mathlab

##### 1 Comment

Ahmed J. Abougarair
on 24 Mar 2024

clc;

clear all;

F = @(t,y) 4*exp(0.8*t)-0.5*y

t0=input('Enter the value of t0 : \n');

y0=input('Enter the value of y0 : \n');

tn=input('Enter the value of t for which you want to find the value of y : \n');

h=input('Enter the step length : \n');

i=0;

while i<tn

k_1 = F(t0,y0);

k_2 = F(t0+0.5*h,y0+0.5*h*k_1);

k_3 = F((t0+0.5*h),(y0+0.5*h*k_2));

k_4 = F(((t0)+h),(y0+k_3*h));

nexty = y0 + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;

y0=nexty;

t0=t0+h;

i=i+h;

end

fprintf('The value of y at t=%f is %f \n',t0,y0)

% validate using a decent ODE integrator

tspan = [0,1]; Y0 = 2;

[tx,yx] = ode45(F, tspan, Y0);

fprintf('The true value of y at t=%f is %f \n',tspan(end),yx(end))

Et= (abs(yx(end)-y0)/yx(end))*100;

fprintf('The value of error Et at t=%f is %f%% \n',tspan(end),Et)

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!