## Search in the document preview

**EXPERIMENT 7
**

**SOLVING ORDINARY DIFFERENTIAL EQUATION
**

**1.0 OBJECTIVES
**

1.. To understand fundamentals of ordinary differential equations 2.. To apply various method for solving ordinary differential equations numerically using

MATLAB software

**2.0. EQUIPMENT
**

Computers and Matlab program in the Mechanical Design.

**3.0 INTRODUCTION & THEORY
An equation involving the derivatives or differentials of the dependent
**

**variable is called a differential equation. A differential equation involving only one
independent is called an ordinary differential equation. If a differential equation
involves two or more independent variables, it is called a partial differential
equation.
**

**Ordinary differential equations are classified according to their order,
linearity, and boundary conditions. The order of an ordinary differential equation is
defined to be the order of the highest derivative present in that equation. Some
examples of first-, second-, and third-order differential equations are
**

**- 1st-order equation
- 2nd-order equation
**

**and
- 3rd-order equation
**

**where x is the independent variable; y is the dependent variable. Ordinary
differential equations can be classified as linear and nonlinear equations. A
differential equation is linear if it can be written in form
**

**Ordinary differential equations can be classified as initial value problems or
boundary value problems. An equation is called an initial value problem (IVP) if the
values of the dependent variables or derivatives are known at the initial value of the
independent variables. An equation for which the values of the dependent variable
or their derivatives are known at the final value of the independent variable is called
a final value problem. If the dependent variable or its derivatives are known at more
than one point of the independent variable, the differential equation is a boundary-
value problem.
**

NT 258 Numerical Analysis Laborat

PAGE 25

**Figure 1: The sequence of events in the application of ODEs for
engineering problem solving
**

**3.1. Ordinary differential equations with MATLAB function
The standard MATLAB package has built-in functions for solving ODEs. The
standard ODE solvers include two functions to implement the adaptive step size
Runge-Kutta method. These are ode23, which uses second- and third-order formula
to attain medium accuracy, and ode45, which uses forth- and fifth-order formulas to
attain higher accuracy.
The full syntax of these functions:
**

*[t, y]=solver (‘your_function’, tspan, y0)
*y Dependent variable

**t Independent variable
**

**Solver ode23 or ode45
**

**your_function This function contains the ODEs you want to solve
**

**tspan A vector specifying the interval of integration, [to, tf]. To obtain
**

**solutions at specific times (all increasing or all decreasing), use
**

**tspan=[to,t1,..,tf].
**

**y0 A vector of initial conditions
**

**There are three ways to define a function
**

**• M-file
**

**• Inline function
**

**Example 7.1
Find the solution of the problem
**

** and exact solution is
in interval 0 ≤ t ≤ 2 with step size is 0.2, using MATLAB functions ode23 and ode45.
**

NT 258 Numerical Analysis Laborat

PAGE 25

**Solution
To define an inline function, by typing:
>> f=inline('t./y')
The other way is via a function M-file, by typing
function z = func(t,y)
z = t./y
To calculate and plot the approximate solution on the interval [0, 2]. Type
>> ode45(f,[0 2],1) - if f is an inline function
>>ode45(@func,[0 2],1) – if func.m is an M-file
There is still another possibility, namely to put the definition of the anonymous
function directly into ode45 command like this:
>>ode45(@(t,y) t./y,[0 2],1)
Plot a family of approximate solutions by using a vector of initial values as the third
argument to ode45 where type
>>ode45(f,[0 2],1:0.2:3)
To extract numerical values type (0de45 does not produce a graph)
>>[t,ya]=ode45(f,[0 2],1)
>>plot(t,ya) – produces a graph without circles
To specify the t values, use a vector (with more than two elements) as the second
argument to ode45 where t = 0, 0.2, 0.4,…,2, type
>>[t,ya]=ode45(f, 0:0.2:2, 1)
To make a table of the values of t and ya(t) with the corresponding values of the
exact solution , type
>> y=sqrt(t.^2+1);
>> format long
>> [t ya y]
ans =
**

** 0 1.00000000000000 1.00000000000000
0.20000000000000 1.01980390030003 1.01980390271856
0.40000000000000 1.07703295800176 1.07703296142690
0.60000000000000 1.16619037435052 1.16619037896906
0.80000000000000 1.28062484165690 1.28062484748657
1.00000000000000 1.41421355599242 1.41421356237310
1.20000000000000 1.56204992884201 1.56204993518133
1.40000000000000 1.72046504740386 1.72046505340853
1.60000000000000 1.88679622083269 1.88679622641132
1.80000000000000 2.05912602304385 2.05912602819740
2.00000000000000 2.23606797273725 2.23606797749979
**

NT 258 Numerical Analysis Laborat

PAGE 25

**Note: **To use **ode23** just replace ‘**ode45**’ with ‘**ode23**’.

**3.2. Numerical Method
3.2.1. Euler’s Method
**Although Euler’s method can be derived in several ways, the derivation based on Taylor’s
series expansion is considered here. The value of can be expressed using Taylor’s series
expansion about *xi*, as

(7.1)

where the third term on the right-hand side of Eq.(7.1) denotes the error or remainder term, and

By Substituting into Eq.(7.1), Eq.(7.1) can be expressed as

(7.2)

If *h* is small, the error term can be neglected, and Eq.(7.2) yields

(7.3) is known as Euler’s or Euler-Cauchy or the point slope method.

**Euler’s Method Procedure
**1. Start with the known initial condition, .
2. Choose *h *.

3. Calculate the solution at each node point by using Eq. (7.3)

4. Check the stability of the solution by changing* h* (make it smaller)
a) If the solution changes for the selected number of significant figures recalculate by

decreasing *h* (go to step 2)
b) If the solution does not change for the selected number of significant figures, that is

the finals

NT 258 Numerical Analysis Laborat

PAGE 25

**Flowchart Euler
**

NT 258 Numerical Analysis Laborat

PAGE 25

**Example 7.2
**Use Euler’s method to integrate from t=0 to 4 with step size of 1. The initial condition at
t=0 is y=2. Note that the exact solution can be determined analytically as

**Procedures-MATLAB Program
**1. Start a new *MatLab* file by clicking on **File**, **New** and **M-file** that opens an empty file in

the Editor/Debugger window.

2. Write the program given below in **M-file
**

%Input - f is the function entered as a string 'f' % - 'a' and 'b' are the left and right endpoints % - 'ya' is the initial condition y(a) – dependent value % - 'dt’ or ‘dt' is step size % - ‘ta’ is the initial condition t(a) – independent value

%Output - E

**df=input('Enter the ordinary differential equation(df):');
f=input('Enter the exact equation(f):');
a=input('Enter the left endpoints(a):');
b=input('Enter the right endpoints(b):');
dt=input('Enter the step size(dt):');
ya=input('Enter the initial condition(ya):');
xa=input('Enter the initial condition(xa):');
**

**M=(b-a)/dt;
t=zeros(M,1);
y_euler=zeros(M,1);
y_exact=zeros(M,1);
error=zeros(M,1);
y_euler(1)=ya;
y_exact(1)=ya;
t(1)=ta;
**

**for j=1:M
t(j+1)=t(j)+dt;
y_euler(j+1)=y_euler(j)+(dt*feval(df,t(j),y_euler(j)));
y_exact(j+1)=feval(f,t(j+1));
error(j+1)=abs((y_exact(j+1)-y_euler(j+1))/y_excat(j+1))*100;
end
E=[t y_euler y_exact error]
**

3. Click on **Save As** to save it as **Euler.m**.
4. To see how it works, type **Euler **in *MatLab Command Window.
*

**3.2.2. Improvements of Euler’s Method
**

**3.2.2.1 Heun’s Method
**

NT 258 Numerical Analysis Laborat

PAGE 25

In Euler’s method, the value of the function *f* (*x*,*y*), which denotes the derivative , is
computed at the beginning of the interval *h* and is assumed to be a constant over the
entire interval. This assumption is a major source of error since the derivative, , changes
from point over the interval *h*. In Heun’s method, the derivative or slope, , is compute at
two points-one at the beginning and the other at the end of the interval h-and their average
value is used to achieve an improvement.

Recall that in Euler’s method, the slope at the beginning of an interval (7.2)

Is used to extrapolate linearly to : (7.3)

For the standard Euler method we would stop at this point. However, in Heun’s method the
calculated in Eq.(7.3) is not the final answer, but an intermediate prediction. This is why we
have distinguished it with a superscript 0. Equation(7.3) is called a *predictor equation*. It
provides an estimate of that allows the calculation of an estimated slope at the end of the
interval:

(7.4)

Thus, the two slopes [Eqs.(7.2) and (7.4)] can be combined to obtain an average slope for the interval:

This average slope is then used to extrapolate linearly from to using Euler’s method:

which is called a *corrector equation*. The Heun method is a *predictor-corrector approach*.

The computational procedure of Heun’s method can be stated as follows:
1. Start with the known initial condition, .
2. Choose *h
*3. Find the value of as

4. Determine the value of *y* at *xi* + 1
-**without iteration
**

or

**-with iteration
**5. Check the convergence criterion as

-**with iteration
**

NT 258 Numerical Analysis Laborat

PAGE 25

6. Check the stability of the solution by changing *h* (make it smaller)
a) If the solution changes for the selected number of significant figures,

recalculate by decreasing *h* (go to step 2)
b) If the solution does not change for the selected number of significant

figures, that is the final solution

**FlowChart of Heun’s Method
**

NT 258 Numerical Analysis Laborat

PAGE 25

`

**Example 7.3
**Use Huen’s method to integrate from t=0 to 4 with step size of 1. The initial condition at
t=0 is y=2. Note that the exact solution can be determined analytically as

**Procedures-MATLAB Program
**

NT 258 Numerical Analysis Laborat

PAGE 25

.1 Start a new *MatLab* file by clicking on **File**, **New** and **M-file** that opens an empty file in
the Editor/Debugger window.

.2 Write the program given below in **M-file
**

%Input - f is the function entered as a string 'f' % - 'a' and 'b' are the left and right endpoints % - 'ya' is the initial condition y(a) % - 'dt' is the step size % - ‘ta’ is the initial condition t(a) %Output – R

**df=input('Enter the ordinary differential equation(df):');
f=input('Enter the exact equation(f):');
a=input('Enter the left endpoints(a):');
b=input('Enter the right endpoints(b):');
dt=input('Enter the step size(dt):');
ya=input('Enter the initial condition(ya):');
xa=input('Enter the initial condition(xa):');
esp=input(‘Enter the percent tolerance(esp):’);
N=input(‘Enter the number of iteration(N):’);
**

**M=(b-a)/dt;
t=zeros(M,1);
y_huen=zeros(M,1);
y_exact=zeros(M,1);
error=zeros(M,1);
y_huen(1)=ya;
y_exact(1)=ya;
t(1)=ta;
**

**for k=1:M
t(k+1)=t(k)+dt;
y_exact(k+1)=feval(f,t(k+1));
y_huen(k+1)=y_huen(k)+feval(df,t(k),y_huen(k))*dt;
%y_huen(k+1)=y_huen(k)+(dt/2*(feval(df,t(k),y_huen(k))+feval
(df,t(k+1),y_huen(k+1))));
**

**for i=1:N
y_huen(k+1)=y_huen(k)+(dt/2*(feval(df,t(k),y_huen(k))+feval
(df,t(k+1),y_huen(k+1))));
error1(k+1)=abs((y_huen(k+1)-y_huen(k))/y_huen(k+1))*100;
if error1 < esp
break
end
**

** y_huen(k+1)=y_huen(k+1)
**

**end
**

** error(k+1)=abs((y_exact(k+1)-y_huen(k+1))/y_exact(k+1))*100
end
**

NT 258 Numerical Analysis Laborat

PAGE 25

**R=[t y_huen y_exact error];
**

.3 Click on **Save As** to save it as **Huen.m**.

.4 To see how it works, type **Huen **in *MatLab Command Window.
*

**3.2.2. Runge-Kutta’s Method
**

These methods are used for the solution of first order ODE (linear and nonlinear). They are a particular set of self- starting numerical methods. They can be used to generate an entire solution. They are more accurate than Euler’s method, but the calculations are more involved

Ruge-Kutta methods require only one initial point to start the procedure. The solution using Ruge-Kutta’s method can be stated in the form:

(7.5)

where is called the increment function, which is chosen to represent the average slope over the interval . The increment function can be expressed as

where *n* denotes the order of the Runge-Kutta’s method; are constants; and are
recurrence relations given by

**.
**

**.
**

**.
**

and (7.6)

where p and a are constants.

**3.2.2.1. Second-Order Runge-Kutta Method
**

We now consider three of most commonly used versions of second-order Runge-Kutta method – Heun’s method, Midpoint’s method, and Ralston’s method.

**Midpoint’s Method
**

(7.7)

with

NT 258 Numerical Analysis Laborat

PAGE 25

**Heun’s Method
**

(7.8)

with

**Ralston’s Method
**

(7.9)

with

**Runge-Kutta’s Method Procedure
**

1. Start with the known initial condition, .

2. Choose *h
*

3. Calculate the solution at each node point by using Eq. (7.7)-** Midpoint’s Method,
**Eq. (7.8)-** Heun’s Method, **or Eq. (7.9)-** Ralston’s Method.
**

4. Check the stability of the solution by changing *h* (make it smaller)
c) If the solution changes for the selected number of significant figures, recalculate

by decreasing *h* (go to step 2)
d) If the solution does not change for the selected number of significant figures, that is

the final solution

NT 258 Numerical Analysis Laborat

PAGE 25

**FlowChart Second-Order Runge-Kutta (Huen’s Method)
**

NT 258 Numerical Analysis Laborat

PAGE 25

**Example 7.3
**Use Euler’s method to integrate from t=0 to 4 with step size of 1. The initial condition at
t=0 is y=2. Note that the exact solution can be determined analytically as

**Procedures-MATLAB Program
**1. Start a new *MatLab* file by clicking on **File**, **New** and **M-file** that opens an empty file in

the Editor/Debugger window.

2. Write the program given below in **M-file
function R=Rk(f,a,b,ya, xa, h)
**

%Input - f is the function entered as a string 'f' % - 'a' and 'b' are the left and right endpoints % - 'ya' is the initial condition y(a) % - 'h' is the step size % - ‘ta’ is the initial condition t(a) %Output – R

**df=input('Enter the ordinary differential equation(df):');
f=input('Enter the exact equation(f):');
a=input('Enter the left endpoints(a):');
b=input('Enter the right endpoints(b):');
dt=input('Enter the step size(dt):');
ya=input('Enter the initial condition(ya):');
ta= input('Enter the initial condition(ta):');
**

**M=(b-a)/dt;
t=zeros(M,1);
y_rk2=zeros(M,1);
y_exact=zeros(M,1);
error=zeros(M,1);
y_rk2(1)=ya;
y_exact(1)=ya;
t(1)=ta;
**

**for k=1:M
t(k+1)=t(k)+dt;
**

NT 258 Numerical Analysis Laborat

PAGE 25

** y_exact(k+1)=feval(f,t(k+1));
k1=feval(df,t(k),y_rk2(k));
k2=feval(df,t(k)+dt,y_rk2(k)+(k1*dt));
y_rk2(k+1)=y_rk2(k)+((k1+k2)/2);
error(k+1)=abs(y_exact(k+1)-y_rk2(k+1)/y_exact(k+1))*100;
end
**

**R=[t y_rk2 y_exact error];
**

3. Click on **Save As** to save it as **Rk.m**.
4. To see how it works, type **Rk2 **in *MatLab Command Window.
*

**3.2.2.2. Fourth-Order Runge-Kutta Methods
**The fourth-order Runge-Kutta methods are most popularly used and more accurate is
compared with the second-order Runge-Kutta methods. The iterative process is given
by (7.10)

where

and

**Runge-Kutta’s Method Procedure
**

1. Start with the known initial condition, .

2. Choose *h*.

3. Calculate the solution at each node point by using Eq. (7.10)

4. Check the stability of the solution by changing *h* (make it smaller)
a). If the solution changes for the selected number of significant figures, recalculate

by decreasing *h* (go to step 2)
b). If the solution does not change for the selected number of significant figures, that is

the final solution

**Flowchart Fourth-Order Runge-Kutta
**

NT 258 Numerical Analysis Laborat

PAGE 25

**Example 7.4
**Use Euler’s method to integrate from t=0 to 4 with step size of 1. The initial condition at
t=0 is y=2. Note that the exact solution can be determined analytically as

NT 258 Numerical Analysis Laborat

PAGE 25

**Procedures-MATLAB Program
**1 Start a new *MatLab* file by clicking on **File**, **New** and **M-file** that opens an empty file in

the Editor/Debugger window.

2 Write the program given below in **M-file
**%Input - f is the function entered as a string 'f'
% - 'a' and 'b' are the left and right endpoints
% - 'ya' is the initial condition y(a)
% - 'dt' is the step size
% - ‘xa’ is the initial condition x(a)
%Output - R

**a=input('Enter the left endpoints:')
b=input('Enter the right endpoints:')
h=input('Enter the step size:')
ya=input('Enter the initial condition y(a):')
ta=input('Enter the initial condition t(a):')
**

**M=(b-a)/dt;
t=zeros(M,1);
y_4rk=zeros(M,1);
y_exact=zeros(M,1);
error=zeros(M,1);
y_4rk(1)=ya;
y_exact(1)=ya;
t(1)=ta;
**

**for k=1:M
t(k+1)=t(k)+dt;
y_exact(k+1)=feval(f,t(k+1));
k1=feval(df,t(k),y_4rk(k));
k2=feval(df,t(k)+(1/2*dt),y_4rk(k)+(1/2*k1*dt));
k3=feval(df,t(k)+(1/2*dt),y_4rk(k)+(1/2*k2*dt));
k4=feval(df,t(k)+dt,y_4rk(k)+(k3*dt));
y_4rk(k+1)=y_4rk(k)+((dt/6)*(k1+(2*k2)+(2*k3)+k4));
error(k+1)=abs((y_exact(k+1)-y_4rk(k+1))/y_exact(k+1))*100;
end
R=[t y_4rk y_exact error]
**3 Click on **Save As** to save it as **Rk4.m**.
4 To see how it works, type **Rk4** in MATLAB**Command Window**.

**3.2.2.3. Higher-Order Runge-Kutta Methods
** (7.10)

NT 258 Numerical Analysis Laborat

PAGE 25

where

and

**Higher-Order Runge-Kutta Methods**
1. Start with the known initial condition, .

2. Choose *h*.
3. Calculate the solution at each node point by using Eq. (7.10).4. Check the stability of the solution by changing* h* (make it smaller)

e) If the solution changes for the selected number of significant figures, recalculate
by decreasing *h* (go to step 2)

f) If the solution does not change for the selected number of significant figures, that is the final solution.

**Flowchart Higher-Order Runge-Kutta Methods
**

NT 258 Numerical Analysis Laborat

PAGE 25

**Example 7.5
**Use Euler’s method to integrate from t=0 to 4 with step size of 1. The initial condition at
t=0 is y=2. Note that the exact solution can be determined analytically as

**Procedures-MATLAB Program
**1. Start a new *MatLab* file by clicking on **File**, **New** and **M-file** that opens an empty file

in the Editor/Debugger window.

2. Write the program given below in

%Input - f is the function entered as a string 'f' % - 'a' and 'b' are the left and right endpoints % - 'ya' is the initial condition y(a) % - 'dt' is the step size % - ‘ta’ is the initial condition t(a) %Output - R

**a=input('Enter the left endpoints:')
b=input('Enter the right endpoints:')
h=input('Enter the step size:')
ya=input('Enter the initial condition y(a):')
ta=input('Enter the initial condition t(a):')
**

**M=(b-a)/dt;
t=zeros(M,1);
y_rf=zeros(M,1);
y_exact=zeros(M,1);
**

NT 258 Numerical Analysis Laborat

PAGE 25

**error=zeros(M,1);
y_rf(1)=ya;
y_exact(1)=ya;
t(1)=ta;
**

**for k=1:M
t(k+1)=t(k)+dt;
y_exact(k+1)=feval(f,t(k+1));
%y_exact(k+1)=2*exp(x(k+1).^2/2)-1;
k1=feval(df,t(k),y_rf(k));
k2=feval(df,t(k)+(1/4*dt),y_rf(k)+((1/4)*k1*dt));
k3=feval(df,t(k)+(1/4*dt),y_rf(k)+(1/8*k1*dt)+(1/8*k2*dt));
k4=feval(df,t(k)+(1/2*dt),y_rf(k)-(1/2*k2*dt)+(k3*dt));
k5=feval(df,t(k)+(3/4*dt),y_rf(k)+(3/16*k1*dt)+(9/16*k4*dt));
k6=feval(df,t(k)+dt,y_rf(k)-(3/7*k1*dt)+(2/7*k2*dt)+(12/7*k3*dt)- (12/7*k4*dt)
+(8/7*k5*dt));
y_rf(k+1)=y_rf(k)+((dt/90)*((7*k1)+(32*k3)+(12*k4)+(32*k5)+(7*k6)));
error(k+1)=abs((y_exact(k+1)-y_rf(k+1))/y_exact(k+1))*100
end
F=[t y_rf y_exact error];
**

3. Click on **Save As** to save it as **Rk.m**.
4. To see how it works, type **Rk** in *MATLAB Command Window.
*

NT 258 Numerical Analysis Laborat

PAGE 25

**4.0 LAB ASSIGNMENT
Assignment
**The free-fall velocity of a parachutist can be simulated as

where *v *= velocity(m/s), *t* = time (s), *g* = acceleration due to gravity (9.81 m/s2), *c*d = drag
coefficient (kg/m), and m = mass (kg). For a 80-kg parachutist, solve this equation using

a) Euler’s Method b) Huen’s Method(Improvement of Euler’s Method) c) Higher Order Runge Kutta’s Method

from t = 0 to 30 s given that v(0) = 0 with step size, h = 2. During free fall, *c*d = 0.25 kg/m.

**Data Analysis
**1. Fill in Table 1, 2 & 3.
2. Plot a graph for results comparison in Figure 1.

**5.0 DATA & RESULTS
**

**Lab Name :** …………………..
**Pc Number :** …………………..
**Folder Name :** …………………..

**Mathematical Model
**

NT 258 Numerical Analysis Laborat

PAGE 25

**Table 1** - **Euler’s Method
**

**K(iteration) v (initial guest) v (Euler)
m/s
**

**v (exact)
m/s
**

**True Percent
Relative Error
**

**(%)
**1
.
.
.
.
.
.
.
.
.
.
n

**Table 2** - **Huen’s Method
**

**K(iteration) v (initial guest) v (Huen)
m/s
**

**v (exact)
m/s
**

**True Percent
Relative Error
**

**(%)
**1
.
.
.

NT 258 Numerical Analysis Laborat

PAGE 25

.

.

.

.

.

.

. n

**Table 3** - **Higher Order Runge Kutta’s Method
**

**K(iteration) v (initial guest) v (Higher Order)
m/s
**

**v (exact)
m/s
**

**True Percent
Relative Error
**

**(%)
**1
.
.
.
.
.
.
.
.
.
.
n

NT 258 Numerical Analysis Laborat

PAGE 25

**6.0 DISCUSSION
**

**7.0 CONCLUSION
**

NT 258 Numerical Analysis Laborat

PAGE 25