Runge-Kutta 4th order method

% 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?

[Naveen]

what can be done if functions are very large?

[Torsten]

Define them as a function instead of a function handle:

https://uk.mathworks.com/help/matlab/ref/function.html

Answer 1

[David Wilson]

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, '--')

[Rashmi Radha]

dear David Wilson

what is use of Y in this code

[Yixuan Qi]

Hi I got "Unable to perform assignment because the left and right sides have a different number of elements." for the y(i+1) line. Any suggestion on how to fix it?

[Walter Roberson]

the code would need to be adjusted slightly if the ode function has more than one state (and so returns a vector.)

[Presley]

Thanks, it was helpful

Answer 2

[Sandip Das]

%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);

Answer 3

[mahmoud mohamed abd el kader]

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

[RITIK PANKAJ]

how can we enter the Input like what would be the format of input

[soham roy]

What modifications do we need to make in this code to solve 3 ODEs with different initial conditions?

[Walter Roberson]

@soham roy

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,:)

Answer 4

[Mj]

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=

[Mj]

????

Answer 5

[David Wilson]

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_a(t) = 

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

f2_b(t) = 

Da = diff(a);
D2a = diff(a,2);

% Now attempt to solve analytically 
dsolve(D2a == f2_b, a(0) == 0.8, Da(0) == 1)

Warning: Unable to find symbolic solution.

 
ans =
 
[ empty sym ]
 

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)

     2
1   a     1
- - -- - ---
2    6   3 a

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)

Answer 6

[Amr Mohamed]

how can we write the code for this problem :

[Moneeb Ur Rehman]

get the y on other side, integrate then to find 1st derivative. Now apply R.k method to solve. Hope you understood;

[Amr Mohamed]

Thanks sir

Answer 7

[monsef]

y=x^2-2yx

h=0.2

y0=0

x0=1

wriet program im mathlab

[Ahmed J. Abougarair]

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)