MATLAB Answers

Runge-Kutta 4th order method

5,137 views (last 30 days)
Show older comments
Mariam Gasra
Mariam Gasra on 5 May 2019
  • 1
Commented: Rashmi Radhaon 17 Feb 2022
% 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
fori=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?

Sign in to answer this question.

Accepted Answer

David Wilson
David Wilson on 6 May 2019
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.5;% 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
fori=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 = -0.5;
[tx, yx] = ode45(F_xy, tspan, y0)
plot(x,y,'o-', tx, yx,'--')
1 Comment
Rashmi Radha
Rashmi Radha on 17 Feb 2022
dear David Wilson
what is use of Y in this code

Sign in to comment.

More Answers (6)

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 ;
fori = 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
1 Comment
RITIK PANKAJ
RITIK PANKAJ on 12 Dec 2020
how can we enter the Input like what would be the format of input

Sign in to comment.

Mj
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=
1 Comment
Mj
Mj on 9 Nov 2020
????

Sign in to comment.

Meysam Mahooti
Meysam Mahooti on 5 May 2021
//www.tatmou.com/matlabcentral/fileexchange/55430-runge-kutta-4th-order?s_tid=srchtitle
David Wilson
David Wilson on 6 May 2021
Ran in:
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.
symsFu c2real
symsa(t)
f2 = (((((2 + c2) * (Fu ^ 2)) / (1 + c2)) + 1) * (^ 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:
symsFu c2real
symsa real
f2 = (((((2 + c2) * (Fu ^ 2)) / (1 + c2)) + 1) * (^ 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)
Amr Mohamed
Amr Mohamed on 9 May 2021
how can we write the code for this problem :
2 Comments
Show 1 older comment
Amr Mohamed
Amr Mohamed on 8 Jun 2021
Thanks sir

Sign in to comment.

Sandip Das
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
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);

Categories

Tags

Community Treasure Hunt

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

Start Hunting!