Transcript Slide 1

Scientific Computing
Numerical Solution Of
Ordinary Differential Equations
- Euler’s Method
Bungee Jumper – Coupled ODEs
Bungee Jumper – Coupled ODEs
• We want to model the vertical dynamics of a jumper connected
to a stationary platform with a bungee cord. F=ma
• Forces:
– mg (gravity, g = acceleration due to gravity)
– cd v2 (drag force, cd = drag coefficient, v = velocity)
(need to always retard v, so if falling (v>0) need force neg,
if rising (v<0) need force pos to reduce dv/dt)
(use sign(v) for drag force)
– k (x-L) (spring force, x = distance measured down from
platform, L = rest length of cord)
– γv
(damping force, γ is damping coefficient of cord)
Bungee Jumper – Coupled ODEs
• We get
dv
k =  = 0 if
2
m
 mg sign(v)cd v  k(x  L)  v x  L
dt
• Dividing by m we get:
cd 2 k
dv

 g  sign( v ) v  ( x  L)  v
dt
m
m
m
k =  = 0 if
xL
Bungee Jumper – Coupled ODEs
• Need to solve two simultaneous ODEs for x and v
 dx
 dt  v

 dv  g  sign( v ) c d v 2  k ( x  L)   v
 dt
m
m
m
k =  = 0 if
xL
Solving First Order ODE’s
• A general form for a first order ODE is
 dx 
F t, x,  0
 dt 
• Or alternatively dx
dt
 f t, x 
 • We desire a solution x(t) which satisfies this
ODE and one specified boundary condition.
x(a) = c

Initial Value Problem
• An Initial Value Problem (IVP) is to solve a
differential equation (of some order) subject
to one or more initial conditions.
Analytic Solutions
• Suppose dx
dt
t
2
x(0) 1
• This can be clearly solved by integration:

3
t
2
 dx   t dt  x  3  C
x(0)  1  C  1
t3
x  1
3
Numerical Solution
• For many ODE’s integration is not possible,
need to do a numerical approximation to the
solution.
• To find a numerical solution we divide the
interval [a,b] for the independent variable t
into n subintervals or steps.


Numerical Solution
• Then the value of the true solution is
approximated at these n values of t.
t i  t 0  ih
i  0,1,...,n
ba
h
n
• We denote the approximation at these pts by
so x i that
x i  xti 
Numerical Solution
• Let x(t) be the actual solution for dx/dt at the
step values.
• Then, assuming no roundoff error, the
difference in the calculated and true value is
the truncation error,
i  xi  xti 

Numerical Solution
Numerical algorithms for solving 1st order odes
for an initial condition are based on one of
two approaches:
1. Direct or indirect use of the Taylor series
expansion for the solution function xt
2. Use of open or closed integration formulas.

Numerical Solution
The various solution algorithms are classified
into:
• One-step: calculation of xi+1 given the
differential equation and previous value: xi
•
Multi-step: in addition to the previous
information the algorithm requires values of
t and x outside of the interval under
consideration.
Numerical Solution
• The first method considered is based on the
Taylor series expansion.
• It forms the basis for some of the other
methods.
Taylor Series
• For this method we express the solution x(t)
about some starting point t0 using a Taylor
expansion.
h2
h3
x t 0  h   x t 0  hx'(t 0 ) 
x''(t 0 ) 
x'''(t 0 )  ...
2
6
• The ODE is given by dx
 f t, x 
dt
Taylor Series
• Thus, by substitution from


x t 0  h   x t 0  hf t 0 , x t 0 
dx
 f t, x 
dt
h2
h3

f  t 0 , x
t0  
f  t 0 , x t 0   ...

2
6




d
f t, xt 
• Where f t, xt  denotes
2 dt
d
and
f '' t, x t  denotes
2 f t, x t 
dt
… etc

Error Analysis for Taylor
Approximation
• Suppose we use the nth order Taylor series
approximation
x(t1 )  xt0  h   xt0   hf t0 , xt0 
h2
h n ( n)

f t0 , xt0    
f t0 , xt0 
2
n!
to generate the value of x1=x(t1). How close are
we to the actual value?
Error Analysis for Taylor
Approximation
• In general, suppose we generate points for x
by:
x (ti 1 )  x ti  h   x ti   hf ti , x ti 
h2
h n (n)

f ti , x ti    
f ti , x ti 
2
n!
to generate the value of xi+1=x(ti+1). How close
are we to the actual value?
Error Analysis for Taylor
Approximation
• The exact value for x(ti+1) is given by


x t i1   x t i  h   x t i  hf t i , x t i 





 ...
h 2 f  t i , x t i 
h 3 f  t i , x t i 
2!
3!

 h
h n f n 1 t i , x t i 
n!
n 1
 in xi , xi 1 
f n  , x 
n 1!
Error Analysis for Taylor
Approximation
• Algorithms for which the last term in
expansion is dropped are of order hn.
• The error is of order hn+1.
• The local truncation error is bounded as
n 1
h
follows et 
M
n 1!
• where

int i ,t i1 
M  f  n  , x 
max



Taylor Series Solution
• Note: differentiation of f(t,x) can be
complicated.
• Direct Taylor expansion is not used other than
the simplest case,


xti1  xt i  hf t i, x t i   Oh
2

Euler’s Method


xti1  xt i  hf t i, x t i 
Euler’s Method
• The problem with the simple Euler method is
the inherent inaccuracy in the formula.
Euler’s Method
• The geometric interpretation is shown the
diagram.
x
x1
x(t)
x(t0)
x(t1)
h
t0
t1
t
Euler’s Method
• The solution across the interval [t0,t1] is
assumed to follow the line tangent to x(t) at
t0.
x
x1
x(t)
x(t0)
x(t1)
h
t0
t1
t
Euler’s Method
• When the method is applied repeatedly across
several intervals in sequence, the numerical
solution traces out a polygon segment with
sides of slope fi, i=0,1,2,…,(n-1).
Euler’s Method
• The simple Euler method is a linear
approximation, and is a perfect solution only if
the function is linear (or at least linear in the
interval).
• This is inherently inaccurate.
• Because of this inaccuracy small step sizes are
required when using the algorithm.
Example
A ball at 1200K is allowed to cool down in air at an ambient temperature
of 300K. Assuming heat is lost only due to radiation, the differential
equation for the temperature of the ball is given by
d
 2.2067 10 12  4  81 10 8 ,  0   1200 K
dt

Find the temperature at

t  480 seconds using Euler’s method. Assume a step size of
h  240 seconds.
http://numericalmethods.eng.usf.edu
Solution
Step 1:
d
 2.2067  10 12  4  81  10 8
dt



f t,   2.20671012 4  81108

 i 1   i  f ti ,  i h
1   0  f t0 ,  0 h
 1200 f 0,1200240



 1200  2.2067 1012 12004  81 108 240
1
 1200  4.5579240
 106.09K
is the approximate temperature at t  t1  t0  h  0  240  240
 240  1  106.09K
http://numericalmethods.eng.usf.edu
Solution Cont
Step 2:
For i  1,
t1  240, 1  106.09
 2  1  f t1 ,1 h
 106.09  f 240,106.09240



 106.09   2.20671012 106.094  81108 240
 106.09  0.017595240
 110.32K
 2 is the approximate temperature at t  t2  t1  h  240 240  480
 480   2  110.32K
http://numericalmethods.eng.usf.edu
Solution Cont
The exact solution of the ordinary differential equation is given by the solution
of a non-linear equation as
0.92593 ln
  300
 1.8519 tan 1 0.00333   0.22067 10 3 t  2.9282
  300
The solution to this nonlinear equation at t=480 seconds is
 (480)  647.57K
http://numericalmethods.eng.usf.edu
Comparison of Exact and Numerical
Solutions
Temperature, θ(K)
1400
1200
1000
Exact Solution
800
600
400
h=240
200
0
0
100
200
300
400
500
Time, t(sec)
Figure 3. Comparing exact and Euler’s method
http://numericalmethods.eng.usf.edu
Effect of step size
Table 1. Temperature at 480 seconds as a function of step size, h
Step, h
(480)
Et
|єt|%
480
240
120
60
30
−987.81
110.32
546.77
614.97
632.77
1635.4
537.26
100.80
32.607
14.806
252.54
82.964
15.566
5.0352
2.2864
 (480)  647.57K
(exact)
http://numericalmethods.eng.usf.edu
Comparison with exact results
Temperature, θ(K)
1500
1000
Exact solution
500
h=120
h=240
0
0
-500
100
200
Tim e, t (sec)
300
400
500
h=480
-1000
-1500
Figure 4. Comparison of Euler’s method with exact solution for different step sizes
http://numericalmethods.eng.usf.edu