Download Solving Ordinary Differential Equations: Numerical Methods and MATLAB Applications and more Exercises Applied Mathematics in PDF only on Docsity!
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
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.
Analysis Laborat
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 ode
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
Example 7.
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.
Analysis Laborat
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 x (^) i , as
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
If h is small, the error term can be neglected, and Eq.(7.2) yields
is known as Euler’s or Euler-Cauchy or the point slope method.
Euler’s Method Procedure
- Start with the known initial condition,.
- Choose h.
- Calculate the solution at each node point by using Eq. (7.3)
- 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
Analysis Laborat
Flowchart Euler
Analysis Laborat
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:
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
-with iteration
5. Check the convergence criterion as
Analysis Laborat
- 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
Analysis Laborat
.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)) end*
Analysis Laborat
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:
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
with
Analysis Laborat
FlowChart Second-Order Runge-Kutta (Huen’s Method)
Analysis Laborat
Example 7. 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;
Analysis Laborat
Example 7. 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
Analysis Laborat
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/2dt),y_4rk(k)+(1/2k1*dt));
k3=feval(df,t(k)+(1/2dt),y_4rk(k)+(1/2k2*dt));
k4=feval(df,t(k)+dt,y_4rk(k)+(k3*dt));
y_4rk(k+1)=y_4rk(k)+((dt/6)(k1+(2k2)+(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)
Analysis Laborat
Example 7. 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);
Analysis Laborat
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/4dt),y_rf(k)+((1/4)k1*dt));
k3=feval(df,t(k)+(1/4dt),y_rf(k)+(1/8k1dt)+(1/8k2*dt));
k4=feval(df,t(k)+(1/2dt),y_rf(k)-(1/2k2dt)+(k3dt));
k5=feval(df,t(k)+(3/4dt),y_rf(k)+(3/16k1dt)+(9/16k4*dt));
k6=feval(df,t(k)+dt,y_rf(k)-(3/7k1dt)+(2/7k2dt)+(12/7k3dt)- (12/7k4dt)
+(8/7k5dt));
y_rf(k+1)=y_rf(k)+((dt/90)((7k1)+(32k3)+(12k4)+(32k5)+(7k6)));
error(k+1)=abs((y_exact(k+1)-y_rf(k+1))/y_exact(k+1))*
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.
Analysis Laborat
