Computational Methods in Physics PHYS 3437

1
Computational Methods in Physics PHYS 3437

• Dr Rob Thacker
• Dept of Astronomy Physics (MM-301C)
• thacker_at_ap.smu.ca

2
Todays Lecture

• Numerical solutions of ordinary differential
  equations
• Euler method
• Runge-Kutta
• Adaptive step size selection
• You can skip this lecture if you have a
  significant amount of experience with numerical
  ODEs

3
Simplest approach to ODEs

• To begin consider a general first order ordinary
  differential equation
• Using a Taylor expansion of the solution y(x)
• Which suggests we predict forward function values
  from a starting point x0,y0

4
Crude (forward) Euler solver graphically
y2
y(x)F(x)
y1
y0
h
x0 x1 x2
5
Useful notation

• n step number
• Dx width of interval

6
Numerical errors

• Discretization error
• Error resulting from the computation of
  quantities that have been calculated using an
  approximation to the true solution
• We may neglect higher order terms for example
• This error is sometimes called truncation error
• Unavoidable problem in numerical integration work
• Would occur even in the presence of infinite
  precision
• Round-off error
• Result of finite precision arithmetic

7
Discretization round-off errors

• Suppose the exact solution at a given value of xn
  is (y)F(xn), the (accumulated) discretization
  error is defined as
• Caused by
• approximate formula to estimate yn1
• Input data at start of step do not necessarily
  correspond to exact soln
• Accumulated round-off error is defined as
• Where Yn is the value we actually calculate after
  round-off rather than the true value yn
• So the absolute value of the total error is given
  by the inequality

8
Problems with the forward Euler method

• This method relies upon the derivative at the
  beginning of each interval to predict forward to
  the end of each interval
• Any errors in the solution tend to get amplified
  quickly
• Calculated solution quickly diverges away from
  the true solution as the error grows
• We can precisely calculate the error in a single
  step as follows

9
Local discretization error in Euler Method
The easiest way to calculate the local error is
to use a Taylor expansion


These three terms correspond to Eulers method -
The local discretization error in the
approximation at one step is thus given by
10
Error over a series of intervals
Be VERY careful about the distinction between The
local discretization error ei and the global
error En accumulated over a series of intervals.
In many books/lectures a clear definition of what
error is being discussed isnt always
given. Roughly speaking, if the local error is
O(hn) then the global error will be O(hn-1).
This happens because n(xn-x0)/h and you sum over
n intervals. Unfortunately, calculating En
usually requires fairly complicated derivations
that we dont have time to do.
y0
x0 x1 x2 xn
11
Modified Euler Method (extrapolate-integrate)

• The Euler method is flawed because we use the
  beginning of the interval to predict to the end
• If we could use the midpoint derivative that
  would usually be expected to be a better
  approximation
• If the slope is changing rapidly across an
  interval using the midpoint slope usually gives
  us something closer to the average over the
  interval (of course it doesnt always have to be
  better, just on average)
• However since ymidf(xmid,ymid) and we dont
  know ymid (only y0) we have to estimate it

12
Modified Euler Method

• Let y(x)f(x,y) and y(x0)y0
• Let the step size be given by hx1-x0,
  f0f(x0,y0), x1/2x0h/2, x3/2x1h/2 etc.

13
Modified Euler Method graphically
In this method the slope is estimated at the
mid-point and the y(x) is integrated forward
using that slope.
y(t)
f1/2
y1
ft1/2
y0
Can show that errors in this method are
x0 x1/2 x1
14
Improved Euler Method

• Similar in concept to the trapezoid rule in
  quadrature, by utilizing an average of the start
  and end point derivatives we can better
  approximate the change in y over the interval
• This is the simplest example of so called
  Predictor-Corrector methods, (which includes
  Runge-Kutta methods)

15
Improved Euler method graphically
y(x)
yn
xn xn1
Can show that errors in this method are also
16
Higher order scheme preliminaries

• As a first motivation consider repeating Improved
  E. method but with two substeps

y(x)F(x)
y0
x0 x1/2 x1
Indicates y1 was obtained with two h/2 steps
17
General points about higher order methods

• Higher order accurate schemes will use more
  substep evaluations
• We need to ensure these are chosen optimally so
  as to produce the most accuracy for the least
  number of function evaluations
• Lets go back to the Improved Euler method on a
  series of steps and apply Richardson-Extrapolation
  -esque approach

18
4th order Runge-Kutta

• For a single Improved Euler method we get
• A global evaluation with n steps will produce an
  error En, of order O(h2)
• If we halve the size of h then the global error
  produced by (D) with the smaller step size must
  be ¼ that of the previous evaluation
• That must correspond to the error produced by
  using the formula y1,h/2 (i.e. C) over the n steps

19
Subtraction step defining

• Since we know the global error term produced by
  using y1,h/2 Err(y1,h/2)K(h/2)2Kh2/4 (E)
• Global error term for y1,h Err(y1,h)Kh2 (F)
• Subtract (F) from 4(E) define new y1

This formula will actually be globally 4th order
accurate!
20
Comments

• There are an infinite number of ways of choosing
  slopes and substeps
• Mathematically all these methods are
  parameterized and limited by using appropriate
  Taylor expansions and matching coefficients
• See W. Gear Numerical Initial Value Problems in
  ODEs, 1971, pgs 27-36
• However, one must always make an arbitrary choice
  to derive a specific algorithm
• Many numerical analysts have considered all the
  possible choices and 4th order Runge-Kutta is
  considered a very good method

21
Classical Runge-Kutta method
Equivalent to the five term Taylor expansion
Interpret the sum over interior points as an
average slope. Local discretization error is
proportional to h5, so halving h leads to a 1/32
reduction in the local disc. error. Very widely
used method
No proof of this algebra is extremely lengthy!
22
Runge-Kutta graphically
y0
x0 x1/2
x1
23
Comparison of methods
Lets compare solutions of the following system
x Euler h0.1 Improved Euler h0.1 RK h0.2 RK h0.1 Exact
0 1.0 1.0 1.0 1.0 1.0
0.2 1.5 1.595 2.5016 2.5050062 2.5053299
0.4 4.4774 5.6099494 5.7776358 5.7927853 5.7942260
0.6 8.9038240 12.442193 12.997178 13.047713 13.052522
0.8 17.537495 27.348020 28.980768 29.130609 29.144880
1.0 34.411490 59.938223 64.441579 64.858107 64.897803
24
Convergence with changes in h
Solutions now taken at x1
h Euler Improved Euler RK Exact
0.2 64.441579 64.897803
0.1 34.411490 59.938223 64.858107 64.897803
0.05 45.588400 63.424698 64.894875 64.897803
0.025 53.807866 64.497931 64.897604 64.897803
0.01 60.037126 64.830722 64.897798 64.897803
0.7 error
0.06 error
10-7 error
RK with 10 steps better than improved Euler with
100 steps!
25
(No Transcript)
26
Adapting the step size

• What value of h is optimal?
• If we consider functions that are smooth in some
  regions and rapidly changing in others then h for
  one region is not appropriate in another
• h ought to vary from one region to another
• Lets examine an algorithm to vary h
• Well use the relative error on a given interval
  as the controlling parameter

27
Algorithm to adapt h

• Let efractional error tolerance (say 10-4)
• Let Dxmax-xmin be the domain over which
  yf(x,y) is to be solved
• Initially let h0.1eD
• Use R-K to find y1y(xminh) in one step h
• Next use two R-K steps of size h/2 to find a new
  y1(xminh) (most accurate estimate yet)
• (Now have two values we can apply
  Richardson-Extrapolation to)

28
Algorithm to adapt h

• 6) Compute y1y(xminh) using y1 and y1 using
  R.E.
• Note since error is prop to h4 must use 1/16
  multiplier
• 7) Estimate error
• 8) If err lt 0.1e ? h is too small. Increase by
  1.5 for next step
• If err gt e ? h is too big. Halve h and
  repeat iteration
• If 0.1 e err e ? h is OK.
• 9) When y value is acceptable increment all x
  values to do next zone
• 10) Stop when xxmax

29
Notes

• The R.E. ensures that even though we appear to do
  3 R-K steps per h zone, we actually get a much
  faster convergence rate than R-K with a step of
  h/3
• The adaptive grid is a necessary part of the
  error estimation process
• Without local adaption you always have to choose
  the h for the hardest regions in your problem
• You may well be wasting time in regions that are
  easy to solve
• It doesnt make any sense to use R-K without
  adaptive stepping!

30
Summary

• As with numerical integration, low order methods
  do not exhibit high accuracy
• RK methods provide a good compromise between
  computational cost and stability
• Adaptive RK methods are extremely powerful and
  using convergence tests that use R.E. ensures
  optimal usage of compute cycles

31
Next lecture

• Issues with the adaptive step size algorithm
• Second order ODEs
• Nth order ODEs