Title: Computational Methods in Physics PHYS 3437
1Computational Methods in Physics PHYS 3437
- Dr Rob Thacker
- Dept of Astronomy Physics (MM-301C)
- thacker_at_ap.smu.ca
2Todays 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
3Simplest 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
5Useful notation
- n step number
- Dx width of interval
-
-
-
-
-
-
6Numerical 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
7Discretization 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
8Problems 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
9Local 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
10Error 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
11Modified 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
12Modified 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.
13Modified 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
14Improved 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)
15Improved Euler method graphically
y(x)
yn
xn xn1
Can show that errors in this method are also
16Higher 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
17General 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
184th 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
19Subtraction 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!
20Comments
- 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
21Classical 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!
22Runge-Kutta graphically
y0
x0 x1/2
x1
23Comparison 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
24Convergence 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)
26Adapting 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
27Algorithm 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)
28Algorithm 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
29Notes
- 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!
30Summary
- 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
31Next lecture
- Issues with the adaptive step size algorithm
- Second order ODEs
- Nth order ODEs