E7 Ordinary Differential Equation New - 0

ENT 258 Numerical Analysis
Laboratory Module
EXPERIMENT 7
SOLVING ORDINARY DIFFERENTIAL EQUATION
1.0 OBJECTIVES
1.1. To understand fundamentals of ordinary differential equations
1.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
dy
ay f (x) - 1st-order equation
dx
x2
d 2 y dy
2 y sin x - 2nd-order equation

dx 2 dx
and
d3y
d2y
dy
a
by
c sin x - 3rd-order equation
3
2
dx
dx
dx
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
a n ( x)
dny
d n 1 y
dy
a
(
x
)
... a1 ( x)
a0 ( x) y f ( x)
n 1
n
n 1
dx
dx
dx
Laboratory Module
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.
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 RungeKutta 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)

Dependent variable
Independent variable
Solver
ode23 or ode45
your_function
This function contains the ODEs you want to solve
tspan
Laboratory Module
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
dy
t
; y (0) 1 and exact solution is y (t )
dt
y
t 2 1
in interval 0 t 2 with step size is 0.2, using MATLAB functions ode23 and ode45.
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 y a (t ) 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 y (0) 1,1.2,1.4,....,3, type
>>ode45(f,[0 2],1:0.2:3)
Laboratory Module
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 y (t ) t 2 1 , type
>> y=sqrt(t.^2+1);
>> format long
>> [t ya y]
ans =
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
Note: To use ode23 just replace ode45 with ode23.
Laboratory Module
3.2. Numerical Method

3.2.1. Eulers Method
Although Eulers method can be derived in several ways, the derivation based on Taylors
series expansion is considered here. The value of y i 1 can be expressed using Taylors
series expansion about xi, as
y i 1 y i hy ' i
h 2 ''
y (ci ); xi ci xi 1 ,
2
(7.1)
where the third term on the right-hand side of Eq.(7.1) denotes the error or remainder
dy
'
( xi ).
term, y i y ( x i ), y i 1 y ( x i 1 ), and y i
dx
By Substituting
y ' i f ( xi , y i )
into Eq.(7.1), Eq.(7.1) can be expressed as
y i 1 y i hf ( xi , y i )
h 2 ''
y (ci )
2
(7.2)
If h is small, the error term can be neglected, and Eq.(7.2) yields

y i 1 y i hf ( xi , y i ); i 0,1,2,....
is known as Eulers or Euler-Cauchy or the point slope method.
(7.3)
Eulers Method Procedure

1. Start with the known initial condition, y x 0 y 0 .
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
Laboratory Module
Flowchart Euler
Start
x(1)=x0, y(1)=y0
steps size, h
No
Yes
STOP
i=1
Stabilit
y
i>N
i =i 1= +N
yi 1 y i hf ( xi , y i )
xi 1 xi h
Laboratory Module
Example 7.2
Use Eulers method to integrate y ' 4e 0.8t 0.5 y 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
4
(e 0.8t e 0.5t ) 2e 0.5t
1.3
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]
Laboratory Module
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 Eulers Method
3.2.2.1 Heuns Method
dy
, is
dx
computed at the beginning of the interval h and is assumed to be a constant over the
dy
entire interval. This assumption is a major source of error since the derivative,
,
dx
dy
changes from point over the interval h. In Heuns method, the derivative or slope,
, is
dx
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.
In Eulers method, the value of the function f (x,y), which denotes the derivative
Recall that in Eulers method, the slope at the beginning of an interval
y i f xi , y i
'
(7.2)
Is used to extrapolate linearly to y i 1 :
y i01 y i f xi , y i h
(7.3)
For the standard Euler method we would stop at this point. However, in Heuns method the
y i01 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 y i 1 that allows the calculation of an estimated slope
at the end of the interval:
y i' 1 f xi 1 , y i01
(7.4)
Thus, the two slopes [Eqs.(7.2) and (7.4)] can be combined to obtain an average slope for
the interval:
y'
y i' y i' 1 f xi , y i f xi 1 , y i01
2
2
This average slope is then used to extrapolate linearly from y i to y i 1 using Eulers
method:
y i 1
f xi , y i f xi 1 , y i01
yi
h
2
which is called a corrector equation. The Heun method is a predictor-corrector approach.
Laboratory Module
The computational procedure of Heuns method can be stated as follows:

1.
Start with the known initial condition, y x 0 y 0 .
2.
Choose h
3.
Find the value of y i 1
as
y i 1 yi hf xi , y i
0
4.
Determine the value of y at xi + 1
y i 1
f xi , yi f xi 1 , y i 1 0
yi h
2
-without iteration
or
f xi , y i f xi 1 , yi k1
2
yik11 yi h
5.
; k 0,1,2,..... -with iteration
Check the convergence criterion as
y ik11 y ik1
s -with iteration
y ik11
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
Laboratory Module
FlowChart of Heuns Method
Start
x(1)=x0, y(1)=y0
Number of steps N
No
Yes
STOP
i=0
Stabilit
y
i>M
i = 1+ i
yi 1 yi hf xi , yi
k=0
k=N
k = 1+ k
f xi , y i f xi 1 , y i k1
y ik11 y i h
10
Laboratory Module
y i k11 y i k1
s
k 1
y i 1
No
Yes
xi 1 xi h
Example 7.3
Use Huens method to integrate y ' 4e 0.8t 0.5 y from t=0 to 4 with step size of 1. The
4
(e 0.8t e 0.5t ) 2e 0.5t
1.3
1. Start a new MatLab file by clicking on File, New and M-file that opens an empty file in

%
%
- 'ya' is the initial condition y(a)
%
- 'dt' is the step size
%
- ta is the initial condition t(a)
%Output R
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);
11
Laboratory Module
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
end
error(k+1)=abs((y_exact(k+1)-y_huen(k+1))/y_exact(k+1))*100
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-Kuttas 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 Eulers method, but the
calculations are more involved
Ruge-Kutta methods require only one initial point to start the procedure. The solution
using Ruge-Kuttas method can be stated in the form:
y i 1 y i h ( xi , y i , h)
(7.5)
where ( xi , y i , h) is called the increment function, which is chosen to represent the

average slope over the interval xi x xi 1 . The increment function can be expressed as
( xi , y i , h) c1 k1 c 2 k 2 ... c n k n
12
Laboratory Module
where n denotes the order of the Runge-Kuttas method; c1 , c 2 ,...c n are constants; and
k1 , k 2 ,..., k n are recurrence relations given by
k1 f ( xi , y i ),
k 2 f ( xi p 2 h, y i a 21 k1 ),
k 3 f ( xi p3 h, y i a31 hk1 a32 hk 2 ),
.
.
.
and
k n f ( xi p n h, y i a n1 hk1 a n 2 hk 2 ... a n , n 1 hk n 1 ).
(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
Heuns method, Midpoints method, and Ralstons method.
Midpoints Method
y i 1 y i hk 2
(7.7)
with
k1 f ( x i , y i )
k 2 f ( xi
h
h
, y i k1 )
2
2
Heuns Method
y i 1 y i
1
h ( k1 k 2 )
2
(7.8)
with
k1 f ( x i , y i )
k2 f ( xi h, yi hk1 )
Ralstons Method
y i 1 y i
1
h ( k1 2 k 2 )
3
(7.9)
13
Laboratory Module
with
k1 f ( x i , y i )
3
3
k2 f ( xi h, yi hk1 )
4
4
Runge-Kuttas Method Procedure

2. Choose h
3. Calculate the solution at each node point by using Eq. (7.7)- Midpoints
Eq. (7.8)- Heuns Method, or Eq. (7.9)- Ralstons Method.
Method,

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 numberStart
of significant figures, that is
the final solution
x(1)=x0, y(1)=y0
Steps size, h
No
STOP
Yes
i=1
Stabilit
y
i>N
i=N
k1 f ( x i , y i )
k 2 f ( xi h, y i hk1 )
FlowChart Second-Order Runge-Kutta (Huens Method)
1
y i 1 y i h(k1 k 2 )
2
14
xi 1 xi h
Laboratory Module
i = i +1
Example 7.3
15
Laboratory Module
4
(e 0.8t e 0.5t ) 2e 0.5t
1.3
1.
Start a new MatLab file by clicking on File, New and M-file that opens an empty
file in the Editor/Debugger window.
function R=Rk(f,a,b,ya, xa, h)

%
%
%
- 'h' is the step size
%
%Output R
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;
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.
16
Laboratory Module
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
h
y i 1 y i (k1 2k 2 2k 3 k 4 )
(7.10)
6
where
k1 f ( x i , y i )
k 2 f ( xi
1
1
h, y i k1 h)
2
2
k 3 f ( xi
1
1
h, y i k 2 h )
2
2
and
k 4 f ( xi h, y i k 3 h)
Runge-Kuttas Method Procedure

2. Choose h.
3. Calculate the solution at each node point by using Eq. (7.10)
a). If the solution changes for the selected number of significant figures, recalculate
b). If the solution does not change for the selected number of significant figures, that is
the final solution
Flowchart Fourth-Order Runge-Kutta
Start
17
Laboratory Module
x(1)=x0, y(1)=y0
steps size, h
No
Yes
STOP
i=1
i>M
Stabilit
y
i = i +1
k1 f ( x i , y i )
1
1
h, y i hk1 )
2
2
1
1
k 3 f ( xi h, y i hk 2 )
2
2
k 4 f ( xi h, y i hk 3 )
k 2 f ( xi
y i 1 y i
h
( k1 2 k 2 2 k 3 k 4 )
6
xi 1 xi h
Example 7.4
18
Laboratory Module
4
(e 0.8t e 0.5t ) 2e 0.5t
1.3
1 Start a new MatLab file by clicking on File, New and M-file that opens an empty file in
2 Write the program given below in M-file

%
%
%
%
- 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;
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
4
Click on Save As to save it as Rk4.m.

To see how it works, type Rk4 in MATLAB Command Window.
3.2.2.3. Higher-Order Runge-Kutta Methods
19
y i 1 y i
Laboratory Module
h
(7 k1 32k 3 12k 4 32k 5 7 k 6 )
90
(7.10)
where
k1 f ( x i , y i )
k 2 f ( xi
1
1
h, y i hk1 )
4
4
k 3 f ( xi
1
1
1
h, y i hk1 hk 2 )
4
8
8
k 4 f ( xi
1
1
h, y i hk1 hk 3 )
2
2
k 5 f ( xi
3
3
9
h, y i hk1 hk 4 )
4
16
16
and
k 6 f ( xi h, y i
3
2
12
12
8
hk1 hk 2 hk 3 hk 4 hk 5 )
7
7
7
7
8
Higher-Order Runge-Kutta Methods

2. Choose h.
3. Calculate the solution at each node point by using Eq. (7.10).
e) If the solution changes for the selected number of significant figures, recalculate
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
Start
x(1)=x0, y(1)=y0
20
Laboratory Module
step size, h
No
Yes
STOP
i=1
i>M
Stabilit
y
i = i +1
k1 f ( x i , y i )
1
1
h, y i k1 )
4
4
1
1
1
k 3 f ( x i h, y i k 1 h k 2 h )
4
8
8
1
1
k 4 f ( xi h, y i k 2 h k 3 h)
2
2
3
3
9
k 5 f x i h, y i
k1 h
k4h
4
16
16
k 2 f ( xi
3
2
12
12
8
k 6 f x i h , y i k1 h k 2 h
k3h k 4 h k5 h
7
7
7
7
7
y i 1 y i
h
(7 k1 32k 3 12k 4 32k 5 7 k 6 )
90
xi 1 xi h
Example 7.5
4
(e 0.8t e 0.5t ) 2e 0.5t
1.3
21
Laboratory Module
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

%
%
%
%
%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);
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)=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.
22
4.0
Laboratory Module
LAB ASSIGNMENT
Assignment
The free-fall velocity of a parachutist can be simulated as
c
dv
g d v2
dt
m
where v = velocity(m/s), t = time (s), g = acceleration due to gravity (9.81 m/s 2), cd = drag
coefficient (kg/m), and m = mass (kg). For a 80-kg parachutist, solve this equation using
23
a)
b)
c)
Laboratory Module
Eulers Method
Huens Method(Improvement of Eulers Method)
Higher Order Runge Kuttas Method
from t = 0 to 30 s given that v(0) = 0 with step size, h = 2. During free fall, cd = 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
24
Laboratory Module
Table 1 - Eulers Method

K(iteration)
v (initial guest)
v (Euler)
m/s
v (exact)
m/s
True Percent
Relative Error
(%)
v (Huen)
m/s
v (exact)
m/s
True Percent
Relative Error
(%)
1
.
.
.
.
.
.
.
.
.
.
n
Table 2 - Huens Method

K(iteration)
v (initial guest)
1
.
.
.
.
25
Laboratory Module
.
.
.
.
.
.
n
Table 3 - Higher Order Runge Kuttas Method

K(iteration)
v (initial guest)
v (Higher Order)
m/s
v (exact)
m/s
True Percent
Relative Error
(%)
1
.
.
.
.
.
.
.
.
.
.
n
Figure 1: v (m/s) vs time (s)
26
6.0
DISCUSSION
7.0
CONCLUSION
Laboratory Module
27
Laboratory Module
28

E7 Ordinary Differential Equation New - 0

Uploaded by

Copyright:

Available Formats

E7 Ordinary Differential Equation New - 0

Uploaded by

Document Information

Original Description:

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

E7 Ordinary Differential Equation New - 0

Uploaded by

Copyright:

Available Formats

ENT 258 Numerical Analysis

3.0 INTRODUCTION & THEORY

2 y sin x - 2nd-order equation

ENT 258 Numerical Analysis

Ordinary differential equations can be classified as initial value problems or

Figure 1: The sequence of events in the application of ODEs for

3.1. Ordinary differential equations with MATLAB function

[t, y]=solver (your_function, tspan, y0)

This function contains the ODEs you want to solve

ENT 258 Numerical Analysis

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

A vector of initial conditions

There are three ways to define a function

ENT 258 Numerical Analysis

To extract numerical values type (0de45 does not produce a graph)

0.20000000000000 1.01980390030003 1.01980390271856

Note: To use ode23 just replace ode45 with ode23.

ENT 258 Numerical Analysis

3.2. Numerical Method

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

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

Eulers Method Procedure

ENT 258 Numerical Analysis

ENT 258 Numerical Analysis

%Input - f is the function entered as a string 'f'

ENT 258 Numerical Analysis

3. Click on Save As to save it as Euler.m.

Recall that in Eulers method, the slope at the beginning of an interval

Is used to extrapolate linearly to y i 1 :

y i' y i' 1 f xi , y i f xi 1 , y i01

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

ENT 258 Numerical Analysis

The computational procedure of Heuns method can be stated as follows:

Find the value of y i 1

Determine the value of y at xi + 1

; k 0,1,2,..... -with iteration

Check the convergence criterion as

Check the stability of the solution by changing h (make it smaller)

ENT 258 Numerical Analysis

FlowChart of Heuns Method

ENT 258 Numerical Analysis

%Input - f is the function entered as a string 'f'

ENT 258 Numerical Analysis

R=[t y_huen y_exact error];

where ( xi , y i , h) is called the increment function, which is chosen to represent the

ENT 258 Numerical Analysis

where p and a are constants.

3.2.2.1. Second-Order Runge-Kutta Method

ENT 258 Numerical Analysis

Runge-Kuttas Method Procedure

4. Check the stability of the solution by changing h (make it smaller)

ENT 258 Numerical Analysis

ENT 258 Numerical Analysis

function R=Rk(f,a,b,ya, xa, h)

ENT 258 Numerical Analysis

3.2.2.2. Fourth-Order Runge-Kutta Methods

Runge-Kuttas Method Procedure

Flowchart Fourth-Order Runge-Kutta