Numerical Methods for Ordinary Differential Equations

58 downloads 158 Views 1MB Size Report
Numerical Methods for. Ordinary Differential Equations. Branislav K. Nikolić. Department of Physics and Astronomy, University of Delaware, U.S.A..
Numerical Methods for Ordinary Differential Equations Branislav K. Nikolić

Department of Physics and Astronomy, University of Delaware, U.S.A. PHYS 460/660: Computational Methods of Physics http://www.physics.udel.edu/~bnikolic/teaching/phys660/phys660.html

Ordinary Differential Equations  d  d2 dn F  y , y (t ), 2 y (t ),..., n y (t )  = 0 dt dt  dt  Ordinary: only one independent variable Differential: unknown functions enter into the equation through its derivatives Order: highest derivative in F Degree: exponent of the highest derivative 3

d  Example:  2 y (t )  − y (t ) = 0  dt  2

PHYS 460/660: Numerical Methods for ODE

What is Solution of ODE?

y = y(t ) A problem involving ODE is not completely specified by its equation ODE has to be supplemented with boundary conditions:

y y

t

•Initial value problem: is given at some starting value i , and it is desired to find at some final points f or at some discrete list of points (for example, at tabulated intervals).

t

•Two point bondary value problem: Boundary conditions are specified at more than one ; typically some of the conditions will be specified at i and some at f .

t

t

PHYS 460/660: Numerical Methods for ODE

t

What is Numerical Solution to the Initial Value Problem?

dy (t ) = f ( t , y (t ) ) ; y (t0 ) = y0 dt

A numerical solution to this problem generates sequence of values for the independent variable

t1 , t2 ,…, tn

and a corresponding sequence of values of the dependent variable

y1 , y2 ,…, yn

so that each

y n approximates solution at t n :

y (t n ) ≈ y n , PHYS 460/660: Numerical Methods for ODE

n = 0,1, …

Euler Metod All finite difference methods start from the same conceptual idea: Add small increments to your function corresponding to derivatives (right-hand side of the equations) multiplied by the stepsize. Euler method is an implementation of this idea in the simplest and most direct form.

y E u le r

y

y tru e

∆t t PHYS 460/660: Numerical Methods for ODE

Single-Step Forward Propagation

Euler Algorithm for First-Order ODE dy = f (t , y ) dt

∆ y = f (t , y ) ∆ t

initialize

t1 , y1 ≡ y ( t1 )

do while i ≤ n y i +1 = y i + f ( t i , y ) ∆ t t i +1 = t i + ∆ t end do PHYS 460/660: Numerical Methods for ODE

Step Size Effects in Radioactive Decay t

− dNU NU Analytics: =− ⇒ NU = NU (t = 0)e τ dt τ

Numerics (Euler): dNU NU (∆t) = NU (0) + ∆t + O (∆t)2 dt Ni Ni+1 ≈ Ni − ∆t

(

)

τ

Numerical solution will depend on the step size

∆t

PHYS 460/660: Numerical Methods for ODE

Stability of Euler Algorithm Step size if often limited by the stability criterion:

dy = − ay ⇒ y (0) = 1, y = e − at dt After n Euler steps of size ∆ t :

y n +1 = y n − ay n ∆ t ⇒ y n = (1 − a ∆ t ) n Approximate solution will decay monotonically only if

∆ t ≤ ∆ t max

1 ≡ a

∆t

is small enough:

For a single decaying exponential-like solution (i.e. if there is only one first order equation) the existence of a stability criterion is not a problem because ∆ t has to be small for the reasons of accuracy. PHYS 460/660: Numerical Methods for ODE

Accuracy: Discretization and Roundoff Errors Lε Integrate over interval: L = t f − t0 ⇒ Full Error: Ch + h p

Local:

C Number of steps for roundoff error to be comparable with the discretization error: N ≈ L    Lε 

1 p +1

du  = f (un , tn ) dt  ⇒ LEn = yn+1 − un+1 (tn+1 ) N −1 t un (tn ) = yn  f = f (t ) ⇒ y (t ) = ∫ f (τ ) dτ ≈ ∑ hn f (tn ) t n=0 Global: t LEn = hn f (tn ) − ∫ f (τ ) dτ t N

0

n +1

GEn = yn − y (tn )

n

N −1

Method is of order n iff: n +1

LEn = O(h ) ⇔ LEn ≤ Ch h = tn +1 − tn ≡ ∆t PHYS 460/660: Numerical Methods for ODE

n +1

GEn = ∑ hn f (tn ) − ∫

tN

t0

n =0 N −1

GEn = ∑ LEn n =0

f (τ ) dτ

Global Discretization Error Example Suppose we want to find the solution over the interval [ 0,T ] →Divide the interval into n equal steps so that ∆t = T n This

T  − aT y(T ) = e , yn = 1 − a  n 

n

is a measure of the global truncation error, error i.e., the error over a fixed range in t. It is proportional to the first power of the step size and hence the Euler method is a first order method (do not confuse this with the fact that we + …are applying it in this case to a first order equation).

(aT )2 (aT )3 y(T ) = 1 − aT + − +… 2! 3! n(n − 1) (aT )2 n(n − 1)(n − 2) (aT )3 − yn = 1 − aT + 2 n 2! n3 3! 1 (aT )2 3 (aT )2  1  a∆t y(T ) − yn = − +…+ O  2  ∼ aTe− aT n 2! n 3! 2 n  PHYS 460/660: Numerical Methods for ODE

Reducing Higher Order ODE to System of First Order ODE Solve higher order ODEs by splitting them into sets of first order equations: d2y dy + p ( t ) + q (t ) y = g (t ) 2 dt dt  dz = g (t ) − p (t ) z − q (t ) y  dy  dt ⇒  z= dt  dy = z  dt T here is no unique w ay to do this:  dz  dp ( t )  = g ( t ) + − q ( t )  y  dt dy  dt  z= + p (t ) y ⇒  dt  dy = z − p ( t ) y  dt PHYS 460/660: Numerical Methods for ODE

Example: Realistic Motion of Baseball   d r    2 v m 2 = mg − B2 v + S0 v × ω dt v 2

v − ω r  v 

F d ra g

xi +1 = x i + vix ∆ t

ω v + ω r

B2 x v =v − vv i ∆ t m y i +1 = y i + v iy ∆ t x i +1

 initialize t1 , y (t1 ) do while i ≤ n    yi +1 = yi + f (ti , yi )∆t ti +1 = ti + ∆t end do PHYS 460/660: Numerical Methods for ODE



x i

v iy+1 = v iy − g ∆ t z i +1 = z i + v iz ∆ t v

z i +1

S 0 v xω =v − ∆t m z i

More Realistic Modeling of Air Flow

PHYS 460/660: Numerical Methods for ODE

ODE of Linear Harmonic Oscillator 2

d θ g + sin θ = 0 2 dt l for small θ ⇒ sin θ ≈ θ 2

d θ g + θ = 0, Ω = 2 dt l 2

Etotal

g l

1 2  dθ  1 2 = ml   + mglθ must be conserved! 2  dt  2

PHYS 460/660: Numerical Methods for ODE

Euler Method for Linear Harmonic Oscillator Switch to dimensionless quantities: d 2θ + θ = 0 ⇒ θ = θ0 sin(Ωt + φ ) 2 dt 2 1  dθ  1 2 Etotal =   + θ 2  dt  2 Euler discretization algorithm:

1 2 ω n +1 = ω n − θ n ∆t   2 ω θ E = + n +1 )   total 2 ( n +1 θ n +1 = θ n + ω n ∆t  ⇒  t n +1 = t n + ∆t

 

PHYS 460/660: Numerical Methods for ODE

 Etotal = En 1 + ∆t 2 

(

)

Euler Fails on θ (t )

PHYS 460/660: Numerical Methods for ODE

Euler Fails on

ω (t )

PHYS 460/660: Numerical Methods for ODE

Euler Fails on Phase Space Trajectory

PHYS 460/660: Numerical Methods for ODE

Can We Resurrect Euler by Using Smaller Step Size?

PHYS 460/660: Numerical Methods for ODE

Cromer Fixed Euler Method for LHO ωn +1 = ωn − θ n ∆t  ωn → ωn +1 ⇒ θ n +1 = θ n + ωn +1∆t t = t + ∆t  n +1 n Apparently trivial trick, but:

1 2 2 2 3 En+1 = En + ωn −θn ∆t + O ∆t 2 θ = θ0 sin(t − t0 ), ω = θ0 cos(t − t0 ) 

(

ω2 −θ 2 =θ02 cos2(t −t0 ) PHYS 460/660: Numerical Methods for ODE

)

over a period

( )

=0

From Euler to Higher Order Algorithms

yn +1 = yn + f (tn , yn ) tn +1 = tn + h PHYS 460/660: Numerical Methods for ODE

Mean value theorem

vs.

exact

y(t + ∆t ) = y(t ) + dy dt t ∆t m

Midpoint Method: Second Order Runge-Kutta s1 = f (tn , yn ) h h s2 = f (tn + , yn + s1 ) 2 2  

3

yn +1 = yn + hs2 + O ( h ) tn +1 = tn + h PHYS 460/660: Numerical Methods for ODE

Classical Runge-Kutta Fo ur th

s1 = f (tn , yn )

-o rd er

me th

od

h h s2 = f (tn + , yn + s1 ) 2 2 h h s3 = f (tn + , yn + s2 ) 2 2 s4 = f (tn + h, yn + hs3 )   

h 5 yn +1 = yn + ( s1 + 2 s2 + 2 s3 + s4 ) + O ( h ) 6 tn +1 = tn + h PHYS 460/660: Numerical Methods for ODE

Classical Runge-Kutta F90 Subroutine

PHYS 460/660: Numerical Methods for ODE

General Single-Step Methods Each of the k stages of the algorithm computes slope si by evaluating f (t , y ) for a particular value of t and a value of obtained by taking linear combinations of the previous slopes:

y

i −1

si = f (t n + α i h , y n + h ∑ β i , j s j ), i = 1, … , k j =1

The proposed step is also a linear combination of the slopes: k

y n +1 = y n + h ∑ γ i s i

i =1 Error is estimated from yet another linear combination of the slopes: k

e n +1 = h ∑ δ i s i i =1

The parameters are determined by matching terms in the Taylor series expansion of the slopes → the order of the method is the exponent of the smallest power of h that cannot be matched.

In MATLAB ODE numerical routines are named as odennxx, where nn indicates the order and xx is some special feature of the method.

PHYS 460/660: Numerical Methods for ODE

Example: MATLAB ode23 Function (Bogacki and Shampine BS23 Algorithm) s1 = f (tn , yn ) h h s2 = f (tn + , yn + s1 ) 2 2 3 3 s3 = f (tn + h, yn + hs2 ) 4  4

h  yn +1 = yn + (2 s1 + 3s2 + 4 s3 )  9  e tn +1 = tn + h; s4 = f (tn +1 , yn +1 ) 

h (−5s1 +6s2 +8s3 −9s4) n+1 = 72

PHYS 460/660: Numerical Methods for ODE

Beyond Runge-Kutta Methods Runge-Kutta methods propagates a solution over an interval by combining the information from several Euler-style steps (each involving one evaluation of the right-hand side f’s), and then using the information obtained to match Taylor series expansion up to some higher order. Richardson extrapolation method used the powerful idea of extrapolating computed result to the value that would have been obtained if the stepsize had been very much smaller than it actually was. In particular, extrapolation to zero stepsize is the desired goal – implemented by Burlich-Stoer algorithm. Predictor-corrector methods store the solution along the way, and use those results to extrapolate the solution one step advanced; they correct the extrapolation using derivative information at the new point. PHYS 460/660: Numerical Methods for ODE

Stiff Systems of Differential Equations Stiffness arises in systems of ODE where there are two or more very different scales of the independent variable: du dv  u = 2e − t − e −1000t = 998u + 1998v, = −999u − 1999v  u = 2 y − z ⇒ dt dt ⇒  −t −1000 t v = − y + z = − + v e e     u (0) = 1, v(0) = 0 Follow the variation in the solution on the shortest length scale to maintain stability of the integration even though accuracy requirements allow for a much larger step size → use implicit methods: explicit

y′ = −cy, c > 0 ⇒ yn+1 = yn + ∆tyn′ = (1− c∆t ) yn ∆t > 2 c ⇔ yn →∞ as n →∞ implicit

y′ = −cy ⇒ yn+1 = yn + ∆tyn′ +1 ⇒ yn+1 =

yn 1+ c∆t

PHYS 460/660: Numerical Methods for ODE

Solutions to Stiffness Beyond Implicit Euler To improve higher-order (than Euler, which is first-order) methods use:

Generalizations of Runge-Kutta methods → Rosenbrock methods and Kaps-Rentrop methods. Burlich-Stoer algorithm generalized to Bader-Deuflhard semi-implicit extrapolation method. Predictor-corrector methods generalized to Gear backward differentiation method. PHYS 460/660: Numerical Methods for ODE