In a previous post, we compared the results from various 2nd order Runge-Kutta methods to solve a first order ordinary differential equation. In this post, I am posting the matlab program. It is better to download the program as single quotes in the pasted version do not translate properly when pasted into a mfile editor of MATLAB or see the html version for clarity and sample output .
Do your own testing on a different ODE, a different value of step size, a different initial condition, etc. See the inputs section below that is colored in bold brick.
% Simulation : Comparing the Runge Kutta 2nd order method of
% solving ODEs
% Language : Matlab 2007a
% Authors : Autar Kaw
% Last Revised : July 12, 2008
% Abstract: This program compares results from the
% exact solution to 2nd order Runge-Kutta methods
% of Heun’s method, Ralston’s method, Improved Polygon
% method, and directly using the three terms of Taylor series
clc
clear all
clc
clf
disp(‘This program compares results from the’)
disp(‘exact solution to 2nd order Runge-Kutta methods’)
disp(‘of Heuns method, Ralstons method, Improved Polygon’)
disp(‘ method, and directly using the three terms of Taylor series’)
%INPUTS. If you want to experiment these are only things
% you should and can change. Be sure that the ode has an exact
% solution
% Enter the rhs of the ode of form dy/dx=f(x,y)
fcnstr=’sin(5*x)-0.4*y’ ;
% Initial value of x
x0=0;
% Initial value of y
y0=5;
% Final value of y
xf=5.5;
% number of steps to go from x0 to xf.
% This determines step size h=(xf-x0)/n
n=10;
%REST OF PROGRAM
%Converting the input function to that can be used
f=inline(fcnstr) ;
% EXACT SOLUTION
syms x
eqn=[‘Dy=’ fcnstr]
% exact solution of the ode
exact_solution=dsolve(eqn,’y(0)=5′,’x’)
% geting points for plotting the exact solution
xx=x0:(xf-x0)/100:xf;
yy=subs(exact_solution,x,xx);
yexact=subs(exact_solution,x,xf);
plot(xx,yy,’.’)
hold on
% RUNGE-KUTTA METHODS
h=(xf-x0)/n;
% Heun’s method
a1=0.5;
a2=0.5;
p1=1;
q11=1;
xr=zeros(1,n+1);
yr=zeros(1,n+1);
%Initial values of x and y
xr(1)=x0;
yr(1)=y0;
for i=1:1:n
k1=f(xr(i),yr(i));
k2=f(xr(i)+p1*h,yr(i)+q11*k1*h);
yr(i+1)=yr(i)+(a1*k1+a2*k2)*h;
xr(i+1)=xr(i)+h;
end
%Value of y at x=xf
y_heun=yr(n+1);
% Absolute relative true error for value using Heun’s Method
et_heun=abs((y_heun-yexact)/yexact)*100;
hold on
xlabel(‘x’)
ylabel(‘y’)
title_name=[‘Comparing exact and Runge-Kutta methods with h=’ num2str(h)] ;
title(title_name)
plot(xr,yr, ‘color’,’magenta’,’LineWidth’,2)
% Midpoint Method (also called Improved Polygon Method)
a1=0;
a2=1;
p1=1/2;
q11=1/2;
%Initial values of x and y
xr(1)=x0;
yr(1)=y0;
for i=1:1:n
k1=f(xr(i),yr(i));
k2=f(xr(i)+p1*h,yr(i)+q11*k1*h);
yr(i+1)=yr(i)+(a1*k1+a2*k2)*h;
xr(i+1)=xr(i)+h;
end
%Value of y at x=xf
y_improved=yr(n+1);
% Absolute relative true error for value using Improved Polygon Method
et_improved=abs((y_improved-yexact)/yexact)*100;
hold on
plot(xr,yr,’color’,’red’,’LineWidth’,2)
% Ralston’s method
a1=1/3;
a2=2/3;
p1=3/4;
q11=3/4;
xr(1)=x0;
yr(1)=y0;
for i=1:1:n
k1=f(xr(i),yr(i));
k2=f(xr(i)+p1*h,yr(i)+q11*k1*h);
yr(i+1)=yr(i)+(a1*k1+a2*k2)*h;
xr(i+1)=xr(i)+h;
end
%Value of y at x=xf
y_ralston=yr(n+1);
% Absolute relative true error for value using Ralston’s Method
et_ralston=abs((y_ralston-yexact)/yexact)*100;
hold on
plot(xr,yr,’color’,’green’,’LineWidth’,2)
% Using first three terms of the Taylor series
syms x y;
fs=char(fcnstr);
% fsp=calculating f'(x,y) using chain rule
fsp=diff(fs,x)+diff(fs,y)*fs;
%Initial values of x and y
xr(1)=x0;
yr(1)=y0;
for i=1:1:n
k1=subs(fs,{x,y},{xr(i),yr(i)});
kk1=subs(fsp,{x,y},{xr(i),yr(i)});
yr(i+1)=yr(i)+k1*h+1/2*kk1*h^2;
xr(i+1)=xr(i)+h;
end
%Value of y at x=xf
y_taylor=yr(n+1);
% Absolute relative true error for value using Taylor series
et_taylor=abs((y_taylor-yexact)/yexact)*100;
hold on
plot(xr,yr,’color’,’black’,’LineWidth’,2)
hold off
legend(‘exact’,’heun’,’midpoint’,’ralston’,’taylor’,1)
% THE OUTPUT
fprintf(‘\nAt x = %g ‘,xf)
disp(‘ ‘)
disp(‘_________________________________________________________________’)
disp(‘Method Value Absolute Relative True Error’)
disp(‘_________________________________________________________________’)
fprintf(‘\nExact Solution %g’,yexact)
fprintf(‘\nHeuns Method %g %g ‘,y_heun,et_heun)
fprintf(‘\nImproved method %g %g ‘,y_improved,et_improved)
fprintf(‘\nRalston method %g %g ‘,y_ralston,et_ralston)
fprintf(‘\nTaylor method %g %g ‘,y_taylor,et_taylor)
disp( ‘ ‘)
________________________________________________________________________________________________
This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu
Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.
The Numerical Methods Guy 
Thanks for Runge-kutta scheme. This scheme will useful to verify rohedi’s formula of the exact solution for the following nonlinear first order ODE :
dx/dt = A*x^3 + B*x^2 + C*x
Sorry, recently I often visit to eqworld.ipmnet,ru , and on the site there are several rohedi’s formulas. The following I present the post of another rohedi’s formula beside rohedi’s formula for Bernoulli Differential Equation that I introduced previously at comment post your numerical blog for problem of Homocide Victim.
%%%%%%%%%%%%%%%%%
Rohedi formula for 3rd RI Differential Equation
Post by afasmt on 23 Nov 2008, 21:11
Dear Rio and All,
To solve the lowerest order of the RI-DE (shortened from Republic of Indonesia Differential Equation) :
dx/dt = A*x^3 + B*x^2 + C*x , (1)
for arbitrary the initial values of t0 and x0, of course there are several ways. The common way i.e by separating variables x and t recomended by Janesh that is by taking to the form :
∫ dx/(A*x^3 + B*x^2 + C*x) = ∫dt + c, (2)
Hence, we can use an appropriate of integral formulation for solving it.
You’re right Rio, my AF(A) diagram that I used to create the New Planck Formula for black-body radiation can be applied to solve eq.(2). For learning more please search the paper on http://rohedi.com. Here I give you all the exact solution of the AF(A) diagram for eq.(1) in form :
D*ln[((x0/x)^2}*(A*x^2 + B*x + C)/(A*x0^2 + B*x0 + C)] +
2*B*[ arctanh{(2*A*x+B)/D} – arctanh{(2*A*x0+B)/D} ] + 2*D*C*(t-t0) = 0 , (3)
where D = sqrt(B^2 – 4*A*C). Appear that the solution of the RI-DE dx/dt = A*x^3 + B*x^2 + C*x is in implicit form like solution of dx/dt = a + b*x^3 discussed on viewtopic.php?f=4&t=41 of Chini DE of mechanics. Further I call eq.(3) as Rohedi Formula for the lowerest order of RI-DE.
Please verify the above exact solution using the usual way or by comparing to the results of symbolic softwares such as Matematica, MapleV, and/or Matlab.
Best Regards,
Rohedi (http://rohedi.com)
Physics Department, Faculty of Mathematics and Natural Sciences,
Sepuluh Nopember Institute of Technology (ITS) Surabaya, Indonesia.
%%%%%%%%%%%%%%%%%%
But now I need your help Mr. Autar, if the above ODE is changed to the form of :
dv/dt = A*v^3 + B*v^2 + C*v ,
where v is velocity of an object movement, how to find the position of object every time x(t), because we know that x(t) = ∫v(t)dt, while rohedi’s formula for v’s solution in implicit form.
Thank’s for your help.
Astrid Denaya Lesa.
LikeLike
Dear Mr.Autar Kaw,
Thanks for your post about the Runge-kutta scheme. This scheme will useful to verify rohedi’s formula of the exact solution for the following nonlinear first order ODE :
dx/dt = A*x^3 + B*x^2 + C*x
Sorry Mr Autar, recently I often visit to eqworld.ipmnet,ru , and on the site there are several rohedi’s formulas. The following I present the post of another rohedi’s formula beside rohedi’s formula for Bernoulli Differential Equation that I introduced previously at comment post your numerical blog for problem of Homicide Victim.
%%%%%%%%%%%%%%%%%
Rohedi formula for 3rd RI Differential Equation
Post by afasmt on 23 Nov 2008, 21:11
Dear Rio and All,
To solve the lowerest order of the RI-DE (shortened from Republic of Indonesia Differential Equation) :
dx/dt = A*x^3 + B*x^2 + C*x , (1)
for arbitrary the initial values of t0 and x0, of course there are several ways. The common way i.e by separating variables x and t recomended by Janesh that is by taking to the form :
∫ dx/(A*x^3 + B*x^2 + C*x) = ∫dt + c, (2)
Hence, we can use an appropriate of integral formulation for solving it.
You’re right Rio, my AF(A) diagram that I used to create the New Planck Formula for black-body radiation can be applied to solve eq.(2). For learning more please search the paper on http://rohedi.com. Here I give you all the exact solution of the AF(A) diagram for eq.(1) in form :
D*ln[((x0/x)^2}*(A*x^2 + B*x + C)/(A*x0^2 + B*x0 + C)] +
2*B*[ arctanh{(2*A*x+B)/D} – arctanh{(2*A*x0+B)/D} ] + 2*D*C*(t-t0) = 0 , (3)
where D = sqrt(B^2 – 4*A*C). Appear that the solution of the RI-DE dx/dt = A*x^3 + B*x^2 + C*x is in implicit form like solution of dx/dt = a + b*x^3 discussed on viewtopic.php?f=4&t=41 of Chini DE of mechanics. Further I call eq.(3) as Rohedi Formula for the lowerest order of RI-DE.
Please verify the above exact solution using the usual way or by comparing to the results of symbolic softwares such as Matematica, MapleV, and/or Matlab.
Best Regards,
Rohedi (http://rohedi.com)
Physics Department, Faculty of Mathematics and Natural Sciences,
Sepuluh Nopember Institute of Technology (ITS) Surabaya, Indonesia.
%%%%%%%%%%%%%%%%%%
But now I need your help Mr.Autar, if the above ODE is changed to the form of :
dv/dt = A*v^3 + B*v^2 + C*v ,
where v is velocity of an object movement, how to find the position of object every time x(t), because we know that x(t) = ∫v(t)dt, while rohedi’s formula for v’s solution in implicit form.
Thank’s for your help.
Astrid Denaya Lesa.
LikeLike
@Denaya Lesa and all visitor this blog
I believe you have any experience in solving differential equation both analytically or by using numerical scheme. Now, I give you a challenge to solve the following of a couple first order differential equation :
dy/dt=ay^2+bx^2+cy+dx+e
dx/dt=ay^2+bx^2+cy+dx+e
for all of the coefficients a, b, c, d, and e are constant, and also for arbitrary of the initial values of t0, x0, and y0.
LikeLike
Dear Professor Autar Kaw,
Denaya is very glad, coz your blog of numerical methods guy now can be browsed from math blog’s Barack Obama supported, that is from this address
http://blueollie.wordpress.com/2009/02/05/daily-kos-the-president-strikes-back/#comment-31157
Of course, if you have spare time, Denaya invites you to visit the blog and leaving some comments for my posts.
Thank you for your attention,
Denaya Lesa.
LikeLike
Hii. can u plz provide MATlab coding for sovling 1st order usinG IMPROVED RUNGE KUTA METHOD
LikeLike