# Differential Equation, Formulas and forms for Applied Mathematics. Taj International College

DOC (481 KB)
24 pages
15Number of visits
Description
Notes on Ordinary Differential Equations.
20 points
this document
Preview3 pages / 24

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

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

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 MATLABCommand Window.

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

NT 258 Numerical Analysis Laborat

PAGE 25

where

and

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), cd = 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, 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

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