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**.

Denaya

said: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.

Denaya Lesa

said: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.

Nadya Fermega

said:@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.

Denaya Lesa

said: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.

Pingback: 2010 in review | The Numerical Methods Guy