Group 7 - Numerical Solutions of ODE

CSU Vision CSU Mission Core Values CSU IGA
Transforming lives by CSU is committed to transform the lives Productivity Compassion Competent
of people and communities through high Accessibility Accountability
quality instruction and innovative Relevance
Self-disciplined
Excellence
research, development, production and
Republic of the Philippines

Cagayan State University
COLLEGE OF ENGINEERING
Carig Sur, Tuguegarao City
DEPARTMENT OF CHEMICAL ENGINEERING
Advance Engineering Mathematics for ChE

(ChE 56)
Second Semester 2016 – 2017
Course Topic: NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL

EQUATIONS USING MATLAB ®
Course Activity: RESEARCH
Name of Students: BENITO, ANGELICA JOYCE Z.

CABBADU, QUENNIE S.
CALLUENG, JOMHEL B.
CATLI, ANGELICA L.
CAUILAN, KEICEE
Program: BSChE
Year Level: III
Date Submitted: May 12, 2017
Instructor: Engr. CAESAR P. LLAPITAN Rating: ________
Date Checked: ________
1
TABLE OF CONTENTS
I. Introduction -------------------------------------------------------------------------------------------- 1
II. Theoretical Background
A. Euler’s Method ------------------------------------------------------------------------------------ 2
B. Heun’s Method ---------------------------------------------------------------------------------- 2-3
C. Midpoint Method ---------------------------------------------------------------------------------- 4
D. Runge-Kutta Methods -------------------------------------------------------------------------- 4-8
D.1 Second-order Runge-Kutta Method
D.2 Fourth-order Runge-Kutta Method
III. Numerical Analysis ------------------------------------------------------------------------------ 9-58
A. Euler’s Method ------------------------------------------------------------------------------- 9-12
1. Algorithm
2. Flow Chart
3. M-file
B. Heun’s Method ----------------------------------------------------------------------------- 13-16
1. Algorithm
2. Flow Chart
3. M-file
C. Ralston’s Method -------------------------------------------------------------------------- 16-19
1. Algorithm
2. Flow Chart
3. M-file
D. Runge-Kutta Method ---------------------------------------------------------------------- 19-23
1. Algorithm
2. Flow Chart
3. M-file
E. Application ---------------------------------------------------------------------------------- 24-58
IV. Generalization ------------------------------------------------------------------------------------ 59-60
V. References --------------------------------------------------------------------------------------------- 61
2
I. Introduction
The fundamental laws of physics, mechanics, electricity, and thermodynamics are usually
based on empirical observations that explain variations in physical properties and states of
systems. Rather than describing the state of physical systems directly, the laws are usually
couched in terms of spatial and temporal changes. These laws define mechanisms of change
and when combined with continuity laws for energy, mass, or momentum, differential
equations result. Subsequent integration of these differential equations results in
mathematical functions that describe the spatial and temporal state of a system in terms of
energy, mass, or velocity variations and the integration can be implemented analytically with
calculus or numerically with the computer. However, a purely mathematical analysis is not
possible to solve differential equations that are complex and nonlinear. That’s why computer
simulations and numerical approximations are needed. Hence, as in figure 1, ordinary
differential equations are used in producing mathematical models of these physical laws and
numerical methods are used to find the numerical approximations of their solutions. Here,
five numerical methods for solving ordinary differential equations, which involves only one
independent variable, are presented.
Figure 1: The sequence of events in the development and solution of ODEs
II. Theoretical Background
Ordinary differential equations of the form

dy
 f (t , y ) (1)
dt
are solved by numerical methods in the general form
yi1  yi  h (2)
where the slope  is called an increment function. According to this equation, the slope
estimate of  is used to extrapolate from an old value to a new value over a distance h. Such
approaches are called one-step methods because the value of the increment function is based
on information at a single point i. These one-step methods which can be expressed in the
3
general form of Eq. (2), with the only difference being the manner in which the slope is
estimated will be discussed next.
A. EULER’S METHOD
Figure 2: Graphical depiction of Euler’s method
In Euler’s method, the first derivative provides a direct estimate of the slope at t i :
  f (ti , yi ) (3)
Where f (ti , yi ) is the differential equation evaluated at ti and yi . This estimate can be
substituted into Eq. (2):
yi1  yi  f (ti , yi )h (4)
This formula is referred to as Euler’s method (or the Euler-Cauchy or point-slope
method). A new value of y is predicted using the slope (equal to the first derivative at the
original value of t) to extrapolate linearly over the step size h (Fig. 2).
B. HEUN’S METHOD
(a) (b)
Figure 3: Graphical depiction of Heun’s method. (a)Predictor and (b) corrector.
One method to improve the estimate of the slope involves the determination of two
derivatives for the interval—one at the beginning and another at the end. The two
4
derivatives are then averaged to obtain an improved estimate of the slope for the entire
interval. This approach, called Heun’s method, is depicted graphically in Fig. 3.
Recall that in Euler’s method, the slope at the beginning of an interval
yi'  f (ti , yi ) (5)
is used to extrapolate linearly to yi1 :
yi01  yi  f (ti , yi )h (6)

For the standard Euler method we would stop at this point. However, in Heun’s method
the y 0 i 1 calculated in Eq. (6) is not the final answer, but an intermediate prediction. This
is why we have distinguished it with a superscript 0. Equation (6) is called a predictor
equation. It provides an estimate that allows the calculation of a slope at the end of the
interval:
yi' 1  f (ti 1 , yi01 ) (7)
Thus, the two slopes [Eqs. (5) and (7)] can be combined to obtain an average slope for the
interval:
f (ti , yi )  f (ti 1 , yi01 )
y'  (8)
2
This average slope is then used to extrapolate linearly from yi to yi 1 using Euler’s
method:
f (ti , yi )  f (ti 1 , yi01 )
yi 1  yi  h (9)
2
which is called a corrector equation.
The Heun’s method is a predictor-corrector approach. As just derived, it can be expressed
concisely as
Predictor (Fig. 3a):
yi01  yim  f (ti , yi )h (10)
Corrector (Fig. 3b):
f (ti , yim )  f (ti 1 , yij11 )
yij1  yim  h
2 (11)
( for j  1,2..., m)
5
C. MIDPOINT METHOD
(a) (b)
Figure 4: Graphical depiction of midpoint method. (a) Predictor and (b) corrector.
This technique uses Euler’s method to predict a value of y at the midpoint of the interval
(Fig. 4a):
h
yi 1/ 2  yi  f (ti , yi ) (12)
2
Then, this predicted value is used to calculate a slope at the midpoint:
yi' 1/ 2  f (ti 1/ 2 , yi 1/ 2 ) (13)
which is assumed to represent a valid approximation of the average slope for the entire
interval. This slope is then used to extrapolate linearly from ti to t i 1 (Fig. 4b):
yi1  yi  f (ti1/ 2 , yi1/ 2 )h (14)
D. RUNGE-KUTTA METHODS
Runge-Kutta (RK) methods achieve the accuracy of a Taylor series approach without
requiring the calculation of higher derivatives. Many variations exist but all can be cast in
the generalized form of Eq. (4):
yi 1  yi  h (15)
where  is called an increment function, which can be interpreted as a representative

slope over the interval. The increment function can be written in general form as
  a1 k1  a 2 k 2  ...  a n k n (16)
where the a' s are constants and the k' s are

k1  f (ti , yi )
k 2  f (t i  p1h, yi  q11k1h)
k3  f (ti  p2 h, yi  q21k1h  q22k2 h)

.
6
.
.
k n  f (ti  pn1 , yi  qn1 ,1 k1h  qn1, 2 k 2 h  ...  qn1 , n1 k n1h
where the p' s and q' s are constants. Notice that the k' s are recurrence relationships.
That is, k1 appears in the equation for k 2 , which appears in the equation for k 3 , and so
forth. Because each k is a functional evaluation, this recurrence makes RK methods
efficient for computer calculations. Various types of Runge-Kutta methods can be devised
by employing different numbers of terms in the increment function as specified by n.
D.1 SECOND-ORDER RUNGE-KUTTA METHODS
The second-order version of Eq. (15) is

yi1  yi  (a1k1  a2 k2 )h (17)
Where
k1  f (ti , yi ) (17a)
k2  f (ti  p1h, yi  q11k1h) (17b)
The values for a1 , a2 , p1 and q11 are evaluated by setting Eq. (17) equal to a second-order
Taylor series. By doing this, three equations can be derived to evaluate the four unknown
constants. The three equations are
a1  a2  1 (18)
a2 p1  1 / 2 (19)
a2 p1  1 / 2 (20)
Because we have three equations with four unknowns, these equations are said to be
underdetermined. We, therefore, must assume a value of one of the unknowns to
determine the other three. Suppose that we specify a value for a2 . Then Eqs. (18) through
(20) can be solved simultaneously for
a1  1  a2 (21)
1
p1  q11  (22)
2 a2
Because we can choose an infinite number of values for a 2 , there are an infinite number
of second-order RK methods. Every version would yield exactly the same results if the
solution to the ODE were quadratic, linear, or a constant. However, they yield different
7
results when (as is typically the case) the solution is more complicated. Three of the most
commonly used and preferred versions are presented next.
i. Heun Method without Iteration ( a2  1 / 2 )
If a2 is assumed to be 1 / 2 , Eqs. (21) and (22) can be solved for a1  1 / 2 and

a1  1 / 2  p1  q11  1 . These parameters, when substituted into Eq. (17), yield
1 1 
yi 1  yi   k1  k 2 h (23)
2 2 
Where
k1  f (ti , yi ) (23a)
k2  f (ti  h, yi  k1h) (23b)

Note that is the slope at the beginning of the interval and is the slope at the end
of the interval. Consequently, this second-order Runge-Kutta method is actually
Heun’s technique without iteration of the corrector.
ii. The Midpoint Method a2  1

If a2 is assumed to be 1, then a1  0, p1  q11  1 / 2 , and Eq. (17) becomes
yi1  yi  k 2 h (24)
Where
k1  f (ti , yi ) (24a)
k2  f (ti  h / 2, yi  k1h / 2) (24b)

This is the midpoint method.
iii. Ralston’s Method a2  2 / 3

Ralston (1962) and Ralston and Rabinowitz (1978) determined that choosing
a2  2 / 3 provides a minimum bound on the truncation error for the second-order
RK algorithms. For this version, a1  1 / 3 and p1  q11  3 / 4 , and Eq. (17)
becomes
1 2 
yi 1  yi   k1  k2 h (25)
3 3 
8
Where
k1  f (ti , yi ) (25a)
 3 3 
k 2  f  ti  h, yi  k1h  (25b)
 4 4 
D.2 FOURTH-ORDER RUNGE-KUTTA METHOD
Figure 5: Graphical depiction of the slope estimates comprising the fourth-order RK

method.
The most popular RK methods are fourth order. As with the second-order approaches,
there are an infinite number of versions. The following is the most commonly used form,
and we therefore call it the classical fourth-order RK method:
1
yi 1  yi  (k1  2k2  2k3  k4 )h (26)
6
Where
k1  f (ti , yi ) (26a)
 1 1 
k 2  f  ti  h, yi  k1h  (26b)
 2 2 
 1 1 
k 3  f  t i  h, yi  k 2 h  (26c)
 2 2 
k4  f (ti  h, yi  k3h) (26d)
Notice that for ODEs that are a function of t alone, the classical fourth-order RK method
is similar to Simpson’s 1/3 rule. In addition, the fourth-order RK method is similar to the
Heun approach in that multiple estimates of the slope are developed to come up with an
9
improved average slope for the interval. As depicted in Fig. 5, each of the k’s represents a
slope. Equation (26) then represents a weighted average of these to arrive at the improved
slope. It is certainly possible to develop fifth- and higher-order RK methods. For
example, Butcher’s (1964) fifth-order RK method is written as
1
yi 1  yi  (7k1  32k3  12k 4  32k5  7k 6 )h (27)
90
Where
k1  f (ti , yi ) (27a)
 1 1 
k 2  f  ti  h, yi  k1h  (27b)
 4 4 
 1 1 1 
k3  f  ti  h, yi  k1h  k 2 h  (27c)
 4 8 8 
 1 1 
k 4  f  ti  h, yi  k 2 h  k3 h  (27d)
 2 2 
 3 3 9 
k5  f  (ti  h, yi  k1h  k 4 h  (27e)
 4 16 16 
 3 2 12 12 8 
k6  f  ti  h, yi  k1h  k 2 h  k3 h  k 4 h  k5 h  (27f)
 7 7 7 7 7 
Although the fifth-order version provides more accuracy, notice that six function
evaluations are required. Recall that up through the fourth-order versions, n function
evaluations are required for an nth-order RK method. Interestingly, for orders higher than
four, one or two additional function evaluations are necessary. Because the function
evaluations account for the most computation time, methods of order five and higher are
usually considered relatively less efficient than the fourth-order versions. This is one of
the main reasons for the popularity of the fourth-order RK method.
10
III. Numerical Analysis
Numerical methods are used to find numerical approximations to the solutions of ordinary
differential equations (ODE). Numerical methods discussed above are used to solve the given
ODE. A general algorithm, flow chart and M-file is shown per method and is discussed
briefly to show how the numerical solution is generated using the MATLAB.
A. Euler’s Method
1. General Algorithm
Step 1: Start.
Step 2: Define the conditions.
Step 3: Define the function.
Step 4: Input the values of the variables.
Step 5: Check if the initial interval is greater than or equal to the final interval.
a. If yes, calculate the equation.
b. If no, that is the final solution.
Step 6: Input the equation of the Euler’s method.
yi 1  yi  hf(ti ,yi )
Step 7: Display the output.
Step 8: End.
11
2. General Flow Chart t0= initial time
START y0= corresponding value
of y at t0
tf= upper boundary of time
t0= 0; y0= 0 t (end time)
h= 0; n= 0 n= number of steps to take
tf= 0 y= 0
dy
 f(t, y)
dx
dy
 (1  4 x) y
dx
dydt= (1+4t)
ta(1)=t0
ya(1)=y0
ta(i+1)=ta(i)+h
ta(1)= ta(1)
ya(1)= ya(1)
tf= tf
Is No
ta>=tf?
Yes
yi1  yi  hf(ti ,yi )
Print
ya]]and ta
END
This flow chart explains the series of steps taken by the software MATLAB. The first step
is to start the program. Then, define the conditions such as t0 as the initial time, y0 as the
corresponding values of y at t0, tf as the upper boundary of time t and n as the numbers of
12
steps to be taken. The function must also be defined, clearly showing the independent and
dependent variable. Next, the value defined beforehand must be specified. After that, a
condition is set, letting the software to solve the equation for Euler’s method for computing
the slope given by y i1  yi  hf (t i , yi ) until it will meet the condition ta>=tf and print the
value of dependent variable y for each value of the independent variable t.
3. General M-file
Two general M-files are provided to show two different outputs. One shows the
computation of the values of y using the Euler’s method per increment of the step size
and the other shows only the values of t and y in a form of a table.
a. Eulers.m (title of the M-file) - shows the computation of the values of y using the
Euler’s method per increment of the step size.
% (___________) //set the file name//
% Solve dy/dt = f(t,y)%

function dydt=f(t,y)
dydt='(_____________)' ; //Insert equation//
f=inline(dydt,'t','y') ;
% t0, initial time
t0=t0 ; //input your initial condition t0//
% y0, corresponding value of y at t0
y0=y0; //input the value of y0 with the initial condition t0//
% tf, upper boundary of time t (end time).
tf=tf; //input the last value between the interval//

% n, number of steps to take
n=n; //set the number of steps to take//

% Displays title information
disp(sprintf('\n\n_________________')) //set the title information//
h=(tf-t0)/n ;
disp(sprintf(' h = ( tf - t0 ) / n '))
disp(sprintf(' = ( %g - %g ) / %g ',tf,t0,n))
disp(sprintf(' = %g',h))
ta(1)=t0 ;
ya(1)=y0 ;
for i=1:n
13
disp(sprintf('\nStep %g',i))
disp(sprintf('---------------------------------------------------'))
% Adding Step Size

ta(i+1)=ta(i)+h ;
% Using eulers formula

ya(i+1)=ya(i)+ f(ta(i),ya(i))*h
disp(sprintf('1) Apply the eulers method to estimate y%g',i))

disp(sprintf(' y%g = y%g +f(ta(i),ya(i))*h',i,i-1))
disp(sprintf(' = %g + %g * %g * %g',ya(i),f(ta(i),ya(i)),h))
disp(sprintf(' = %g\n',ya(i+1)))
disp(sprintf(' at t%g = %g',i,ta(i+1)))
end
disp(sprintf('\n\n*********************Results***********************'))
--------------------------------------------------------------------------
b. Euler.m (title of the M-file) - shows only the values of t and y in a form of a table.
function euler(func,y0,dt,t0,tf)
% our time interval

t=t0:dt:tf;
% setting initial y value

y(1)=y0;
% loop using euler's method

for i = 1:length(t)-1
y(i+1) = y(i) + dt*(feval(func,t(i),y(i)));
end
% prints column vectors of t and y

t=t';
y=y';
plot(t,y)
xlabel('Time')
ylabel('y')
disp(t)
disp(y)
14
B. Heun’s Method
Step 1: Start.
Step 6: Input the equation of the Heun’s method.
yi 1  yi 
 
f xi , yi   f xi 1 , yi01
h
2
Step 8: End.
15
2. General Flow Chart
t0= initial time
y0= corresponding value
START of y at t0
t (end time)
t0= 0; y0= 0 n= number of steps to take
h= 0; n= 0
tf= 0 y= 0
ta(1)=t0
ya(1)=y0
ta(i+1)=ta(i)+h
ta(1)= ta(1)
ya(1)= ya(1)
tf=tf
Is No
ta>=tf?
Yes
Print
ya and ta
END
16
condition is set, letting the software to solve the equation for Heuns method for computing
the slope given by yi 1  yi 

 
f xi , yi   f xi 1 , yi01
h until it will meet the condition ta>=tf
2
and print the value of dependent variable y for each value of the independent variable t.
3. General M-file
Heuns.m (title of the M-file)
% (___________) //set the file name//

% t0, initial time


disp(sprintf('\n\n_____________________'))//set the title information//
h=(tf-t0)/n ;
ta(1)=t0 ;
ya(1)=y0 ;
for i=1:n
disp(sprintf('---------------------------------------------------'))
% Adding Step Size

ta(i+1)=ta(i)+h ;
17
% Using heuns formula
ya(i+1)=ya(i)+0.5*h*(f(ta(i),ya(i))+f(ta(i)+h,ya(i)+h*f(ta(i),ya(i))))
disp(sprintf('1) Apply the heuns method to estimate y%g',i))

disp(sprintf(' y%g = y%g +
0.5*h*(f(ta(i),ya(i))+f(ta(i)+h*ya(i)+h*f(ta(i),ya(i))))',i,i-1))
disp(sprintf(' = %g + %g * %g *
%g',ya(i),0.5,h,(f(ta(i),ya(i))+f(ta(i)+h,ya(i)+h*f(ta(i),ya(i)))),h))
end
disp(sprintf('\n\n***********************Results**************************'
))
-------------------------------------------------------------------------
C. Ralston’s Method
Step 1: Start.
a. If yes, calculate the slope
Step 6: Calculate the k’s, each represents a slope.
Step 7: Input the equation of the Ralston method.
1
y i 1  y i  h(k1  2k 2 )
3
Step 9: End.
18
t0= initial time

of y at t0
tf= 0 y= 0
dy
 f(t, y)
dx
dy
 (1  4 x) y
dx
dydt= (1+4t)
ta(1)=t0
ya(1)=y0
ta(i+1)=ta(i)+h
ta(1)= ta(1)
ya(1)= ya(1)
tf= tf
Is No
ta>=tf?
Yes
k1  f ( xi , y i )
3 3
k2  f ( xi  h, yi  hk1 )
4 4
1
y i 1  y i  h(k1  2k 2 )
3
Print
ya and ta
END
19
condition is set, letting the software to solve the equation for Ralston method for computing
1
the slope given by y i 1  y i  h(k1  2k 2 ) until it will meet the condition ta>=tf
3
and print the value of dependent variable y for each value of the independent variable t.
3. General M-file
Ralston.m (title of the M-file)
% (___________) //set the file name//

% t0, initial time


disp(sprintf('\n\n_________________ ')) //set the title information//
h=(tf-t0)/n ;
ta(1)=t0 ;
ya(1)=y0 ;
for i=1:n
disp(sprintf('---------------------------------------------------'))
20
% Adding Step Size
ta(i+1)=ta(i)+h ;
% Calculating k1, and k2
k1 = f(ta(i),ya(i)) ;
k2 = f(ta(i)+0.75*h,ya(i)+0.75*k1*h ) ;
% Using ralston formula
ya(i+1)=ya(i)+((1/3)*k1+(2/3)*k2)*h ;
disp('1) Find k1 and k2 using the previous step information.')

disp(sprintf(' k1 = f( x%g , y%g )',i-1,i-1))
disp(sprintf(' = f( %g , %g )',ta(i),ya(i)))
disp(sprintf(' = %g\n',k1))
disp(sprintf(' k2 = f( x%g + 0.75*h , y%g + 0.75*ki*h )',i-1,i-1))
disp(sprintf(' = f( %g + 0.75 * %g , %g + 0.75 * %g *
%g)',ta(i),h,ya(i),k1,h))
disp(sprintf(' = f( %g , %g )',ta(i)+0.75*h,ya(i)+0.75*k1*h))
disp(sprintf('2) Apply the ralston formula to estimate y%g',i))

disp(sprintf(' y%g = y%g + k2*h',i,i-1))
disp(sprintf(' = %g + %g * %g * %g',ya(i),((1/3)*k1+(2/3)*k2),h))
end
))
D. Runge-Kutta Method
Step 1: Start.
a. If yes, calculate the slope
Step 7: Input the equation of the Runge-kutta method.
ya(i+1)=ya(i)+1/6*(k1+2*k2+2*k3+k4)*h
Step 9: End
21
t0= initial time
START of y at t0
t (end time)
h= 0; n= 0
tf= 0 y= 0
ta(1)=t0
ya(1)=y0
ta(i+1)=ta(i)+h
ta(1)= ta(1)
ya(1)= ya(1)
tf= tf
No
Is
ta>=tf?
Yes
k1 = f(ta(i),ya(i)))
k2 = f(ta(i)+0.5*h,ya(i)+0.5*k1*h)
k3 = f(ta(i)+0.5*h,ya(i)+0.5*k2*h)
k4 = f(ta(i)+h,ya(i)+k3*h)
Print
ya and ta
END
22
condition is set, letting the software to solve the equation for Runge-Kutta method for
1
computing the slope given by ya (i 1)  ya (i )  (k 1  2k 2  2k 3  k 4 )h until it will meet the
6
condition ta>=tf and print the value of dependent variable y for each value of the independent
variable t.
3. General M-file
Two general M-files are provided to show two different outputs. One shows the
computation of the values of y using the Euler’s method per increment of the step size
and the other shows only the values of t and y in a form of a table.
a. Runge-Kutta.m (title of the M-file) - shows the computation of the values of y using the
Euler’s method per increment of the step size.
% (___________) //set the file name//
% Solve dy/dt = f(t,y)

% t0, initial time


disp(sprintf('\n\n__________________')) //set the title information//
h=(tf-t0)/n ;
23
ta(1)=t0 ;
ya(1)=y0 ;
for i=1:n
disp(sprintf('---------------------------------------------------'))
% Adding Step Size

ta(i+1)=ta(i)+h ;
% Calculating k1, k2, k3, and k4
k1 = f(ta(i),ya(i)) ;
k2 = f(ta(i)+0.5*h,ya(i)+0.5*k1*h) ;
k3 = f(ta(i)+0.5*h,ya(i)+0.5*k2*h) ;
k4 = f(ta(i)+h,ya(i)+k3*h) ;
% Using 4th Order Runge-Kutta formula

ya(i+1)=ya(i)+1/6*(k1+2*k2+2*k3+k4)*h ;

disp(sprintf(' k2 = f( x%g + 0.5 * h , y%g + 0.5 * k1 * h )',i-1,i-
1))

1))
disp(sprintf(' k4 = f( x%g + h, y%g + k3*h)',i-1,i-1))
disp(sprintf(' = f( %g , %g )',ta(i)+h,ya(i)+k3*h))
disp(sprintf('2) Apply the Runge-Kutta 4th Order method to estimate

y%g',i))
disp(sprintf(' y%g = y%g + 1/6*( k1 + 2*k2 + 2*k3 +k4) * h',i,i-1))
%g',ya(i),1/6,(k1+2*k2+2*k3+k4),h))
end
disp(sprintf('\n\n*******************Results**************************')
24
b. RK.m (title of the M-file) - shows only the values of t and y in a form of a table.
function yout = ode4(F,t0,h,tfinal,y0)
%ODE4 Classical Runge-Kutta ODE solver.
% yout = ODE4(F,t0,h,tfinal,y0) uses the classical
% Runge-Kutta method with fixed step size h on the interval
% t0 <= t <= tfinal
% to solve
% dy/dt = F(t,y)
% with y(t0) = y0.
y = y0;
yout = y;
for t = t0 : h : tfinal-h
s1 = F(t,y);
s2 = F(t+h/2,y+h*s1/2);
s3 = F(t+h/2,y+h*s2/2);
s4 = F(t+h,y+h*s3);
y = y + h*(s1 + 2*s2 + 2*s3 + s4)/6;
yout = [yout; y];
end
25
E. Application
PROBLEM 1:
Solve the following problem over the interval form x = 0 to 1 using a step size of 0.25
where y(0)= 1. Display your results on graph.
dy
 (1  4 x) y
dx
a) Analytically.
b) Using Euler’s method.
c) Using Heun’s method without iteration.
d) Using Ralston’s method.
e) Using Midpoint Method.
f) Using the fourth-order RK method.
SOLUTION:
A. Analytical Method
Given the ordinary differential equation,

dy
 (1  4 x) y
dx
The analytical solution of the above equation is evaluated as
dy
1
 (1  4 x)dx
2
y
1
y dy  (1  4 x)dx
2
Integrating
1
1
y2 4x 2
 x C
1 2
1
2
In simplifying we obtain
1
2 y  x  2x 2  C
2
Using the initial condition y (0) =1, we solve for the value of the constant of integration C
26
1
2(1) 2  0  2(0) 2  C
C2
Therefore, the particular solution is
( x  2 x 2  2) 2
y
4
B. Euler’s Method
The Euler’s method is is given as

yi 1  yi  f (t i , yi )h
Applying the Euler’s method for h = 25 with the initial condition, y = 1 at t = 0, the
equation will be
y (1)  y (0)  f (0,1)(1)
The slope estimate at t=0 and y = 1 is
f (0,1)  (1  4(0))( 1)
f (0,1)  1
Therefore,
y (1)  1  1(0.25)
y (1)  1.25
The true solution at t=1 is
( x  2 x 2  2) 2
y
4
(1  2(1) 2  2) 2
y
4
y  6.25
27
Euler’s Method Algorithm
Step 1: Start.
t0= 0; y0= 0
h= 0; n= 0
tf= 0 y= 0
dy
 f(t, y)  (1  4t ) y
dx
ta(1)= 0
ya(1)= 1
tf=1
Step 6: Input the equation of the Euler’s method.
yi 1  yi  hf(ti ,yi )
Step 8: End.
28
Euler’s Method Flowchart
t0= initial time

of y at t0
tf= 0 y= 0
dy
 f(t, y)  (1  4t ) y
dx
dy
 (1  4 x) y
dx
ta(1)=t0
dydt= (1+4t)
ya(1)=y0
ta(i+1)=ta(i)+h
ta(1)= 0
ya(1)= 1
tf=1
Is No
ta>=tf?
Yes
yi1  yi  hf(ti ,yi )
Print
ya]]and ta
END
29
MATLAB Program
Mfile for the solution of ODE using Euler’s Method showing the computation of the values
of per increment of the step size
Eulers.m – Title of the Mfile
% eulers method

dydt='(1+4*t)*sqrt(y)' ;
% t0, initial time
t0=0 ;
y0=1;
tf=1;
n=4;
disp(sprintf('\n\nthe eulers method of solving differential equations'))
h=(tf-t0)/n ;
ta(1)=t0 ;
ya(1)=y0 ;
for i=1:n
disp(sprintf('-----------------------------------------------------'))
% Adding Step Size

ta(i+1)=ta(i)+h ;
% Using eulers formula

ya(i+1)=ya(i)+ f(ta(i),ya(i))*h
disp(sprintf('1) Apply the eulers method to estimate y%g',i))

disp(sprintf(' y%g = y%g +f(ta(i),ya(i))*h',i,i-1))
disp(sprintf(' = %g + %g * %g * %g',ya(i),f(ta(i),ya(i)),h))
end
30
disp(sprintf('\n\n*********************Results***********************')
Output
31
32
Mfile for the solution of ODE using Euler’s method showing only the values of t and y:
Euler.m – Title of the Mfile
function euler(func,y0,dt,t0,tf)
% our time interval

t=t0:dt:tf;
% setting initial y value

y(1)=y0;
% loop using euler's method

for i = 1:length(t)-1
y(i+1) = y(i) + dt*(feval(func,t(i),y(i)));
end
% prints column vectors of t and y

t=t';
y=y';
plot(t,y)
xlabel('Time')
ylabel('y')
33
disp(t)
disp(y)
Output
TABLE 1: Comparison of the numerical solution computed using Euler’s method with the
true solution
x ytrue yeuler %Error
0 1.00 1.00 0
0.25 1.41 1.25 12.8
0.5 2.25 1.80902 24.38
0.75 3.7539 2.81776 33.17
1 6.25 4.49638 39.00
Figure 6: True Solution vs. Numerical Solution using Euler’s Method
34
7
4
y
Analytical
3
Euler's
2
0
0 0.2 0.4 0.6 0.8 1 1.2
x
Shown in table 1 are the true and numerical values of the integral of y'  (1  4 x) y ,
with the initial condition y(0) = 1. The numerical values were computed using Euler’s
method with a step size of 0.25. The true solution at t = 0.25 is 1.41 and the numerical
solution is 1.25, therefore, the true percent relative error is 12.8%. Moreover, shown in figure
6 is the graph of the true and numerical solution. The graph shows that the computation
captures the general trend of the true solution so the error is considerable. However, this error
can be reduced by using a smaller step size.
C. Heun’s Method
First, the slope is at (t0, y0) is
y' 0  (1  4(0))( 1)
y'0  1
Then, the predictor is used to compute the value at 1.0:
y '1  1  1(0.25)
y '1  1.25
35
Now, to improve the estimate the value of y i 1 we use the value y ' 0 to predict the slope at
the end of the interval
y'1  f ( xi , y'0 )  (1  4(1)) 1.25
y'1  5.59017
Which can be combine to the initial slope to yield an average slope over the interval from
t=0 to 1
1.25  5.59017
y' 
2
y '  3.42
This result can be substituted to the corrector to give the prediction at t=0.25:
y '1  1  3.42(0.25)
y '1  1.855
Heun’s Method Algorithm
Step 1: Start.
t0= 0; y0= 0
h= 0; n= 0
tf= 0 y= 0
36
dy
 f(t, y)  (1  4t ) y
dx
ta(1)= 0
ya(1)= 1
tf=1
Step 6: Input the equation of the Heuns method.
f (ti , yi )  f (ti 1 , yi01 )

yi 1  yi  h
2
Step 8: End.
37
Heun’s Method Flowchart
t0= initial time
of y at t0
START
t (end time)
n= number of steps to take
t0= 0; y0= 0
h= 0; n= 0
tf= 0 y= 0
ta(1)=t0
ya(1)=y0
ta(i+1)=ta(i)+h
ta(1)= 0
ya(1)= 1
tf=1
Is
No
ta>=tf?
Yes
Print
ya and ta
END
38
MATLAB Program
Mfile for the soltution of ODE using Heun’s Method:
Heuns.m – title of the M-file
% heuns method

% t0, initial time
t0=0 ;
y0=1;
tf=1;
n=4;
disp(sprintf('\n\nthe heuns method of solving differential equations'))
h=(tf-t0)/n ;
ta(1)=t0 ;
ya(1)=y0 ;
for i=1:n
disp(sprintf('---------------------------------------------------'))
% Adding Step Size

ta(i+1)=ta(i)+h ;
% Using heuns formula

ya(i+1)=ya(i)+0.5*h*(f(ta(i),ya(i))+f(ta(i)+h,ya(i)+h*f(ta(i),ya(i))))
disp(sprintf('1) Apply the heuns method to estimate y%g',i))

disp(sprintf(' y%g = y%g +
0.5*h*(f(ta(i),ya(i))+f(ta(i)+h*ya(i)+h*f(ta(i),ya(i))))',i,i-1))
%g',ya(i),0.5,h,(f(ta(i),ya(i))+f(ta(i)+h,ya(i)+h*f(ta(i),ya(i)))),h))
end
39
))
Output
40
41
TABLE 2: Comparison of the numerical solution computed using Heun’s method with the
true solution
x ytrue yheuns %Error
0 1.00 1.00 0
0.25 1.41 1.40451 0.39
0.5 2.25 2.23073 0.86
0.75 3.7539 3.70609 1.24
1 6.25 6.15178 1.6
Figure 7: True Solution vs. Numerical Solution using Heun’s Method

7
4
y
3 Analytical
Heun's
2
0
0 0.2 0.4 0.6 0.8 1 1.2
x
with the initial condition y(0) = 1, where the numerical values are computed using Heun’s
solution is 1.4045, therefore, the true percent relative error is 0.39% and the error is lesser
than the error in the Euler’s method. Moreover, shown in figure 7 is the graph of the true and
numerical solution. The graph shows that the computation captures the trend of the true
solution and the error is very small. Decreasing the step size decreases the error at a faster
rate than for Euler’s method. Thus, the Heun’s method is more accurate than the Euler’s
method.
42
D. Ralston’s Method
k1  f (0,1)  (1  4(0)) 1
k1  1
k 2  f (t i  p1h, yi  q11k1 h)
k 2  f (0.1875,1.875)  1  4(0.1875)( 1.1875 )
k 2  1.9070
1 2 
y1  1   (1)  (1.9070)  x0.25
3 3 
y1  1.40117
Ralston Method Algorithm
Step 1: Start.
t0= 0; y0= 0
h= 0; n= 0
tf= 0 y= 0
dy
 f(t, y)  (1  4t ) y
dx
ta(1)= 0
ya(1)= 1
tf=1
Step 7: Input the equation of the Ralston’s method.
1
y i 1  y i  h( k1  2k 2 )
3
Step 9: End.
43
Ralston Method Flowchart
t0= initial time

of y at t0
tf= 0 y= 0
dy
 (1  4 x) y
dx
dydt= (1+4t) y
ta(1)=t0
ya(1)=y0
ta(i+1)=ta(i)+h
ta(1)= 0
ya(1)= 1
tf=1
Is No
ta>=tf?
Yes
k1  f ( xi , y i )
3 3
k2  f ( xi  h, yi  hk1 )
4 4
1
y i 1  y i  h(k1  2k 2 )
3
Print
ya and ta
END
44
MATLAB Program
Mfile for the soltution of ODE using Midpoint Method:
Ralston.m – title of the M-file
% ralston method

% t0, initial time
t0=0 ;
y0=1;
tf=1;
n=4;
disp(sprintf('\n\nthe ralston method of solving differential equations'))
h=(tf-t0)/n ;
ta(1)=t0 ;
ya(1)=y0 ;
for i=1:n
disp(sprintf('-----------------------------------------------------'))
% Adding Step Size

ta(i+1)=ta(i)+h ;
% Calculating k1, and k2
k1 = f(ta(i),ya(i)) ;
k2 = f(ta(i)+0.75*h,ya(i)+0.75*k1*h ) ;
% Using ralston formula

ya(i+1)=ya(i)+((1/3)*k1+(2/3)*k2)*h ;

disp(sprintf(' k2 = f( x%g + 0.75*h , y%g + 0.75*ki*h )',i-1,i-1))
45
disp(sprintf('2) Apply the ralston formula to estimate y%g',i))

disp(sprintf(' y%g = y%g + k2*h',i,i-1))
disp(sprintf(' = %g + %g * %g * %g',ya(i),((1/3)*k1+(2/3)*k2),h))
end
disp(sprintf('\n\n**********************Results*************************'))
Output
46
47
Table 3: Comparison of the numerical solution computed using Ralston’s method with the
true solution
x ytrue yralston %Error
0 1.00 1.00 0
0.25 1.41 1.40117 0.63
0.5 2.25 2.22012 1.34
0.75 3.7539 3.68678 1.82
1 6.25 6.11935 2.135
Figure 8: True Solution vs. Numerical Solution using Ralston’s Method

7
4
y
3 Analytical
Ralston
2
0
0 0.2 0.4 0.6 0.8 1 1.2
x
48
with the initial condition y(0) = 1, where the numerical values are computed using Ralston’s
solution is 1.4012, therefore, the true percent relative error is 0.63% and the error is lesser
than the error in the Euler’s method but greater than the Heun’s method. Moreover, shown in
figure 8 is the graph of the true and numerical solution. The graph shows that the
computation captures the trend of the true solution and the error is very small. Thus, the
Ralston’s method is more accurate than the Euler’s method but less accurate than the Heun’s
method.
E. Fourth-order RK method
The slope at the beginning of the interval is computed as
k1  f (0,1)  (1  4(0)) 1
k1  1
This value is used to compute a value of y and the slope at the midpoint
y (0.125)  1  1(0.125)
y (0.125)  1.125
k 2  f (0.125,1.125)  (1  4(0.125))( 1.125 )
k 2  1.59099
This slope in turn is used to compute another value y and another slope at the midpoint:
y (0.125)  1  (0.125)(1.59099)
y (0.125)  1.19887
k 3  f (0.125,1.19887)  (1  4(0.125))( 1.19887 )
k 3  1.6423
Next, this slope is used to compute a value of y and a slope at the end of the interval:
y (1)  1  1(0.4106)
y (1)  1.4106
k 4  f (0.25,1.4106)  (1  4(0.25)) 1.4106 )
k 4  2.3753
Finally, the four slope estimates are combined to yield an average slope. This average
slope is then used to make the final prediction at the end of the interval.
49

1
1  2(1.59099)  2(1.6423)  2.3753
6
  1.6403
y (1)  1  1.6403(1)
y (1)  2.6403
Runge-Kutta Method Algorithm
Step 1: Start.
t0= 0; y0= 0
h= 0; n= 0
tf= 0 y= 0
dy
 f(t, y)  (1  4t ) y
dx
ta(1)= 0
ya(1)= 1
tf=1
Step 7: Input the equation of the Runge-Kutta’s method.
1
yi 1  yi  (k1  2k2  2k3  k4 )h
6
Step 9: End.
50
Runge-Kutta Method Flowchart
t0= initial time
START of y at t0
t (end time)
h= 0; n= 0
tf= 0 y= 0
ta(1)=t0
ya(1)=y0
ta(i+1)=ta(i)+h
ta(1)= 0
ya(1)= 1
tf=1
No
Is
ta>=tf?
Yes
k1 = f(ta(i),ya(i)))
k2 = f(ta(i)+0.5*h,ya(i)+0.5*k1*h)
k3 = f(ta(i)+0.5*h,ya(i)+0.5*k2*h)
k4 = f(ta(i)+h,ya(i)+k3*h)
Print
ya and ta
END
51
MATLAB Program
Mfile for the solution of ODE using Runge-Kutta Method showing the computation of the
values of per increment of the step size
Runge-Kutta.m – Title of the Mfile
% RK4 method

% t0, initial time
t0=0 ;
y0=1;
tf=1;
n=4;
disp(sprintf('\n\nThe 4th Order Runge-Kutta Method of Solving Ordinary
Differential Equations'))
h=(tf-t0)/n ;
ta(1)=t0 ;
ya(1)=y0 ;
for i=1:n
disp(sprintf('-----------------------------------------------------------
------'))
% Adding Step Size

ta(i+1)=ta(i)+h ;
% Calculating k1, k2, k3, and k4
k1 = f(ta(i),ya(i)) ;
k2 = f(ta(i)+0.5*h,ya(i)+0.5*k1*h) ;
k3 = f(ta(i)+0.5*h,ya(i)+0.5*k2*h) ;
k4 = f(ta(i)+h,ya(i)+k3*h) ;
% Using 4th Order Runge-Kutta formula

ya(i+1)=ya(i)+1/6*(k1+2*k2+2*k3+k4)*h ;
52
1))

1))
disp(sprintf(' k4 = f( x%g + h, y%g + k3*h)',i-1,i-1))
disp(sprintf(' = f( %g , %g )',ta(i)+h,ya(i)+k3*h))
disp(sprintf('2) Apply the Runge-Kutta 4th Order method to estimate

y%g',i))
disp(sprintf(' y%g = y%g + 1/6*( k1 + 2*k2 + 2*k3 +k4) * h',i,i-1))
%g',ya(i),1/6,(k1+2*k2+2*k3+k4),h))
end
disp(sprintf('\n\n*******************Results**************************')
Output
53
54
55
Mfile for the solution of ODE using Euler’s method showing only the values of t and y
RK.m – Title of the Mfile
function yout = ode4(F,t0,h,tfinal,y0)

%ODE4 Classical Runge-Kutta ODE solver.
% yout = ODE4(F,t0,h,tfinal,y0) uses the classical
% Runge-Kutta method with fixed step size h on the interval
% t0 <= t <= tfinal
% to solve
% dy/dt = F(t,y)
% with y(t0) = y0.
y = y0;
yout = y;
for t = t0 : h : tfinal-h
s1 = F(t,y);
s2 = F(t+h/2,y+h*s1/2);
s3 = F(t+h/2,y+h*s2/2);
s4 = F(t+h,y+h*s3);
y = y + h*(s1 + 2*s2 + 2*s3 + s4)/6;
yout = [yout; y];
end
Output
56
Table 4: Comparison of the numerical solution computed using Runge-Kutta method with the
true solution
x ytrue yrunge-kutta %Error
0 1.00 1.00 0
0.25 1.41 1.41004 0.0028
0.5 2.25 2.24979 0.009
0.75 3.7539 3.75341 0.013
1 6.25 6.24905 0.015
Figure 9: True Solution vs. Numerical Solution using Runge-Kutta Method

7
4
y
3
Analytical
Runge-Kutta
2
0
0 0.2 0.4 0.6 0.8 1 1.2
x
with the initial condition y(0) = 1, where the numerical values are computed using Runge-
Kutta method with a step size of 0.25. The true solution at t = 0.25 is 1.41 and the numerical
solution is 1.41, therefore, the true percent relative error is 0.015%. Of all the four methods,
the Runge-Kutta method gave the least percentage error for the results obtained. Moreover,
shown in figure 9 is the graph of the true and numerical solution. The graph shows that the
computation exactly captures the trend of the true solution and the difference between the
values computed by this method compared to the true values are negligible. Thus, the Runge-
Kutta method is the most accurate among the four numerical methods.
57
PROBLEM 2:
Figure 10: Spherical Tank
A spherical tank has a circular orifice in its bottom through which the liquid flows out
(Fig.10). The flow rate through the hole can be estimated as
Qout  CA 2 gh
where Qout = outflow (m3/s), C = an empirically derived coefficient, A =the area of the orifice
(m2), g =the gravitational constant (=9.81 m/s2), and h= the depth of liquid in the tank. Use
one of the numerical methods described in this chapter to determine how long it will take for
the water to flow out of a 3-m diameter tank with an initial height of 2.75 m. Note that the
orifice has a diameter of 3cm and C=0.55.
SOLUTION:
It is required to determine the time for the water to flow out of the 3-m diameter tank
using the numerical method. The flow rate of the water through the hole is given as
Qout  CA 2 gh
The initial height of the water in the tank is 2.75 m. The above equation can be written as
dV
 CA 2 gh
dt
The volume of the water in the tank is related to the height as
h 3 (3r  h )
V
3
Differentiating the above equation, it is obtained that
dV dh
 ( 2rh  h 2 )
dt dt
58
Substituting the above result in the original equation, it is obtained that
dh
 CA 2 gh  ( 2rh  h 2 )
dt
dh  CA 2 gh

dt 2rh  h 2
The area of the orifice is
A  r 2
A   0.015
2
A  0.000707m 2
Therefore, the required equation to be integrated is
dh 0.550.000707 29.81h

dt 2 1.5h  h 2
dh h
 5.4825  10  4
dt (3h  h 2 )
Integrating
(3h  h 2 )
dh  5.4825  10 4 dt
h
3h h2
 h
dh  
h
dh  5.4825  10 4  dt
5
3
2h 2
2h  2
 5.4825  10 4 t  C
5
Using the initial condition h (0) =2.75 m, we solve for the value of the constant of integration
5
3
2(2.75) 2
2(2.75) 2
 5.4825  10  4 (0)  C
5
C = 4.1043
Therefore, the particular solution is
5
3 2
2h
2h  2
 5.4825  10 4 t  4.1043
5
59
 Using the Runge-Kutta Method
60
Table 5: Numerical Values Computed by the Runge-Kutta Method
t h 4.25 0.2167 8.75 0.1357 13.25 0.1097

0 2.7500 4.50 0.2055 9.00 0.1313 13.50 0.1067
0.25 1.6269 4.75 0.1955 9.25 0.1271 13.75 0.1040
0.50 1.1565 5.00 0.1864 9.50 0.1232 14.0 0.1013
0.75 0.8971 5.25 0.1781 9.75 0.1195 14.25 0.0988
1.00 0.7328 5.50 0.1705 10.0 0.1160 14.50 0.0964
1.25 0.6193 5.75 0.1635 10.25 0.1128 14.75 0.0942
1.50 0.5362 6.00 0.1571 10.50 0.1097 15.0 0.0920
1.75 0.4728 6.25 0.1511 10.75 0.1067 15.25 0.1195
2.00 0.4228 6.50 0.1955 11.0 0.1456 15.50 0.1160
2.25 0.3824 6.75 0.1864 11.25 0.1405 15.75 0.1128
2.50 0.3490 7.00 0.1781 11.50 0.1357 16.0 0.1097
2.75 0.3210 7.25 0.1705 11.75 0.1313 16.25 0.1067
3.00 0.2972 7.50 0.1635 12.0 0.1271 16.50 0.1040
3.25 0.2766 7.75 0.1571 12.25 0.1232 16.75 0.1013
3.50 0.2587 8.00 0.1511 12.50 0.1195 17.0 0.0988
3.75 0.2430 8.25 0.1456 12.75 0.1160 17.25 0.0964
4.00 0.2291 8.50 0.1405 13.0 0.1128 17.5 0.0942
61
IV. Generalization
Table 6: Results of the five different methods used to calculate the solution of the given ODE
with an interval fom x = 0 to 1 using a step size of 0.25
x yanalytical yeuler yheuns yralston yrunge-kutta
0 1.00 1.00 1.00 1.00 1.00
0.25 1.41 1.25 1.40451 1.40117 1.41004
0.5 2.25 1.80902 2.23073 2.22012 2.24979
0.75 3.7539 2.81776 3.70609 3.68678 3.75341
1 6.25 4.49638 6.15178 6.11935 6.24905
Table 7: Percent Relative Error of the Numerical Methods

 euler  heuns  ralston  runge kutta
12.8 0.39 0.63 0.0028
24.38 0.86 1.34 0.009
33.17 1.24 1.82 0.013
39.00 1.6 2.135 0.015
62
Figure 11: Graph of the solution with the given interval from x = 0 to 1 using a step size of
0.25 using five different methods
7
4 Analytical
Euler's
y
3 Heun's
Ralston
2
Runge-Kutta
0
0 0.2 0.4 0.6 0.8 1 1.2
x
The figure shows the values of y for the function y'  (1  4x) y from the interval x=0 to
1with the step size of 0.25. Using five different methods, it can be observed that four out of
the five methods in obtaining the value of y cluster at almost the same points and the other
method came out close hence. Looking at the curves individually, taking the plot of the
function using the analytical method as basis, the curve for Euler’s method, represented with
a red curve, deviate from the basis with the margin of error that ranges from 12.8%-39%. It
can be said observed that the Euler’s method is the least accurate method in solving the
differential equation. Next, the curve for Heun’s method represented with a green curve,
deviate from the basis with the margin of error that ranges from 0.39%-1.6%. The curve
follows the same path with the basis curve and has a very small margin of error. It can be
said that the method is accurate and effective in solving the differential equation. Third, the
curve for the Ralston’s method represented with the violet curve, deviate from the basis curve
with a margin of error that ranges from 0.63%-2.14%. The Ralston’s method can be also
observed as an accurate method in solving the ODE. Lastly, the curve for the Runge-Kutta
represented with a light blue curve, deviate from the basis curve with a margin of error that
ranges from 0.0028%-0.015%. It can be said that the last method is the most accurate method
in solving the ordinary differential equation. It can be observed that as the step size increases,
the value of y obtained increases in percent relative error. Overall, the five different methods
63
showed different margin of error in solving one-step methods in solving ODE. The Runge-
Kutta is taken as the most accurate approach and the Euler’s method as the least accurate
method. The data solved using MATLAB, showed the numerical values needed for the
graphical interpretation. MATLAB is indeed an effective software in solving different
differential equations especially problems related in engineering.
V. References
Chapman, S. J. (2001). MATLAB Programming for Engineers.
Chapra, S. C. (2005). Applied Numerical Methods with MATLAB for Engineers and Scientist.
New York: McGraw-Hill Icn.
Hahn, B. D., & Valentine, D. T. (2001). Essential MATLAB for Engineers and Scientist.
Oxford: Elsvier Ltd.
64

Group 7 - Numerical Solutions of ODE

Uploaded by

Copyright:

Available Formats

Group 7 - Numerical Solutions of ODE

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Group 7 - Numerical Solutions of ODE

Uploaded by

Copyright:

Available Formats

CSU Vision CSU Mission Core Values CSU IGA

Republic of the Philippines

DEPARTMENT OF CHEMICAL ENGINEERING

Advance Engineering Mathematics for ChE

Course Topic: NUMERICAL SOLUTION OF ORDINARY DIFFERENTIAL

Course Activity: RESEARCH

Name of Students: BENITO, ANGELICA JOYCE Z.

Year Level: III

Date Submitted: May 12, 2017

Instructor: Engr. CAESAR P. LLAPITAN Rating: ________

Date Checked: ________

Figure 1: The sequence of events in the development and solution of ODEs

II. Theoretical Background

Ordinary differential equations of the form

Figure 2: Graphical depiction of Euler’s method

Figure 3: Graphical depiction of Heun’s method. (a)Predictor and (b) corrector.

is used to extrapolate linearly to yi1 :

yi01  yi  f (ti , yi )h (6)

yi1  yi  f (ti1/ 2 , yi1/ 2 )h (14)

where  is called an increment function, which can be interpreted as a representative

where the a' s are constants and the k' s are

k3  f (ti  p2 h, yi  q21k1h  q22k2 h)

D.1 SECOND-ORDER RUNGE-KUTTA METHODS

The second-order version of Eq. (15) is

k2  f (ti  p1h, yi  q11k1h) (17b)

i. Heun Method without Iteration ( a2  1 / 2 )

If a2 is assumed to be 1 / 2 , Eqs. (21) and (22) can be solved for a1  1 / 2 and

k2  f (ti  h, yi  k1h) (23b)

ii. The Midpoint Method a2  1

k2  f (ti  h / 2, yi  k1h / 2) (24b)

iii. Ralston’s Method a2  2 / 3

D.2 FOURTH-ORDER RUNGE-KUTTA METHOD

Figure 5: Graphical depiction of the slope estimates comprising the fourth-order RK

yi1  yi  hf(ti ,yi )

% (___________) //set the file name//

% Solve dy/dt = f(t,y)%

% t0, initial time

t0=t0 ; //input your initial condition t0//

% y0, corresponding value of y at t0

y0=y0; //input the value of y0 with the initial condition t0//

% tf, upper boundary of time t (end time).

tf=tf; //input the last value between the interval//

n=n; //set the number of steps to take//

% Adding Step Size

% Using eulers formula

disp(sprintf('1) Apply the eulers method to estimate y%g',i))

% our time interval

% setting initial y value

% loop using euler's method

% prints column vectors of t and y

the slope given by yi 1  yi 

Heuns.m (title of the M-file)

% (___________) //set the file name//

% Solve dy/dt = f(t,y)%

% t0, initial time

t0=t0 ; //input your initial condition t0//

% y0, corresponding value of y at t0

y0=y0; //input the value of y0 with the initial condition t0//

% tf, upper boundary of time t (end time).