Transcript High Performance Computing 811

Computational Methods
in Physics
PHYS 3437
Dr Rob Thacker
Dept of Astronomy & Physics (MM-301C)
[email protected]
Today’s Lecture





Notes on projects & talks
Issues with adaptive step size selection
Second order ODEs
nth order ODEs
Algorithm outline for handling nth order ODEs
Advanced warning – no class on March 31st
Project time-line




By now you should have approximately ½ of the
project coding & development done
There are only 4 weeks until the first talk
I would advise you to have at least some outline
of how your write-up will look already
Project reports are due on the last day of term
Talk Format





Allowing for changing over presenters,
connecting/disconnecting etc, we can allow up
to 60 minutes per lecture period
3 people per lecture means 20 minutes per
presenter
15 minute presentation
5 minutes for questions
Remember your presentation will be marked by
your peers
Issues with the adaptive step size
algorithm

Step (8) If err < 0.1e  h is too small. Increase by 1.5 for next step
If err > e  h is too big. Halve h and repeat iteration
If 0.1 e ≤ err ≤ e  h is OK.


Clearly if there is problems with convergence this step will continue to keep
dividing h forever. You need to set up a limit here, either by not allowing h to
get smaller than some preset limit or counting the number of times you halve
h
Because of rounding error as you increment x it is very unlikely that
you will precisely hit xmax with the final step

Therefore you should choose h=min(h,xmax-x) where x is current position
Issues with the adaptive step size
algorithm: more “gotcha’s”

Step (7) Estimate error
err 


y1*  yˆ1


1 *
y1  yˆ1
2
Should add a very small (“tiny”) number to the absolute value of the
denominator to avoid a divide by zero (when the two estimates are very
close to zero you might - in incredibly rare situations - get the two
numbers adding to zero).
Since you store values of y and x in an array it is possible that
you may run out of storage (because you need too many points).
You should do some check of the number of positions stored
versus the maximum possible allowed in your decleration of the
x & y arrays . This will avoid the possibility of a segmentation
fault.
Second order ODEs

Consider the general second order ODE

We now require two initial values be provided, namely
the y0 value and the derivative y´(x0)
y' ' ( x)  g ( x, y, y' )


These are called the Cauchy conditions
If we let z=y´, then z´=y´´ and we have
y'  z
;
y ( x0 )  y0
z '  y ' '  g ( x, y ( x), z ( x)) ; z ( x0 )  z0
(1)
Second order ODEs cont
  y
  y0 
 z
Let Y    so Y0   , and F   
z
g 
 z0 
then we can write the system as
 
  y0 
 y ' 
 z '   F  Y '  F with Y0   z 
 
 0



We have thus turned a second order ODE into a first order
ODE for a vector
Can apply the R-K solver to the system but you now have two
components to integrate
At each step must update all x, y and z values
Diagrammatically
Remember z’=g(x,y,y’)
Remember y’=z
y(x)
y1
y´0=z0
y0
x0
x1
z´0=g0
z1
z0
x0
x1
z(x)
nth order ODEs



Systems with higher than second order derivatives are
actually quite rare in physics
Nonetheless we can adapt the idea for 2nd order
systems to nth order
Suppose we have a system specified by
y ( x)  g ( x, y, y' , y' ' ,..., y
( n)

( n 1)
)
Such an equation requires n initial values for the
derivatives, suppose again we have the Cauchy conds.
y( x0 )  y0 , y' ( x0 )  y'0 , y' ' ( x0 )  y' '0 ,..., y ( n1) ( x0 )  y ( n1) 0
Definitions for the nth order system
Let
y1  y
y2  y'1
y3  y' 2

yn  y' n-1
y' n  y'' n-1  ...  y2
(n1 )
 y1  y (n)
(n)
Thus y' n  g(x,y,y',y'',...,y (n-1 ) ) and this suggests
writing in a vector format.
Definitions for the nth order system
 y1   y 
 y   y' 
  2  1 
Let Y   y3    y '2  taking derivative s gives

  



  
 yn   y 'n 1 
y2

 y '1  

 y'  
y
3

 2  

   


 
 
 
Y'   y 'i   
yi 1
  F ( x, Y )

   


 

yn

 y 'n 1  
 y '   g ( x, y ,.., y )
n 
1
 n  
So another system of coupled first order equations that can be solved via R-K.
Algorithm for solving such a system
Useful parameters to set:
imax=number of points at which y(x) is to be evaluated
(1000 is reasonable)
nmax=highest order of ODE to be solved (say 9)
errmax=highest tolerated error (say 1.0×10-9)
Declare x(imax),y(imax,nmax),y0(nmax),yl(nmax),
ym(nmax),yr(nmax),ytilde(nmax),ystar(nmax)
Note that the definitions used here are not quite consistent with the
variable definitions used in the discussion of the single function case.
User inputs

Need to get from the user (not necessarily at run
time)
The g(x,y,y’,…) function
 Domain of x, i.e. what are xmin and xmax
 What is the order of the ODE to be solved, stored
in variable nord
 What are the initial values for y and the derivatives
(the Cauchy conditions), these are stored in
y0(nmax)

Correspondence of arrays to variables

The arrays y(imax,nmax), y0(nmax) corresponds to
the y derivatives as follows:



y(1:imax,1)≡y vals
y(1:imax,2)≡y´ vals
y(1:imax,3)≡y´´ vals
|

y(1:imax,nord)≡y(nord-1) vals
with y(1,1)=y0=y0(1)
with y(1,2)=y´0=y0(2)
with y(1,3)=y´´0=y0(3)
|
with y(1,nord)=y(n-1)0 = y0(nord)
Choose initial step size & initialize yl
values

Apply same criterion as standard Runge-Kutta
dx=0.1×errmax×(xmax-xmin)
 dxmin=10-3×dx



x(1)=xl=xmin


We can use this value to ensure that adaptive step size is
never less than 1/1000th of the initial guess
Set initial position for solver
Initialize yl (left y-values for the first interval)

do n=1,nord
yl(n)=y0(n)=y(1,n)
Start adaptive loop

Set i=1 this will count number of x positions evaluated




dxh=dx/2
xr=xl+dx
xm=xl+dxh
--- half width of zone
--- right hand boundary
--- mid point of zone
Perform R-K calculations on this zone for all yn


Need to calculate all y values on right boundary using a single
R-K step (stored in ytilde(nmax))
Need to calculate all y values on right boundary using two
half R-K steps (stored in ystar(nmax)) – for example
call rk(xl,yl,ytilde,nord,dx)
call rk(xl,yl,ym
,nord,dxh)
call rk(xm,ym,ystar ,nord,dxh)
Now evaluate R.E.-esque value yr(n)


err=0.
do n=1,nord
yr (n)  (16. * ystar (n)  ytilde (n)) / 15.
yr (n)  ystar (n)
errn  2
yr (n)  ystar (n)
err  max( err , errn )

Error is now set over all the
functions being integrated
If err< 0.1×errmax then increase dx: dx=dx*1.5
If err > errmax  dx is too big: dx=max(dxh,dxmin) and repeat evaluation
of this zone
If 0.1 errmax ≤ err ≤ errmax  dx is OK & we need to store all results
Update values and prepare for next
step




Increment i: i=i+1 (check that it doesn’t exceed
imax)
x(i)=xr
xl=xr (note should check xl hasn’t gone past
xmax & that h value is chosen appropriately)
do n=1,nord
y(i,n)=yr(n)
 yl(n)=yr(n)


Then return to top of loop and do next step
Runge-Kutta routine



Subroutine rk(xl,yl,yr,nord,dx)
yl and yr will be arrays of size nord
Set useful variables:




dxh=dx/2.
xr=xl+dx
xm=xl+dxh
Recall, given y’=f(x,y), (xl,yl), and dx then
~
~
h
yr  yl  ( f 0  2 f1/ 2  2 f1/ 2  f1 )
6
~
where f 0  f ( xl , yl )
y1/ 2  yl  dxh  f 0
~
~
~
f  f ( x , y ) y  y  dxh  f
1/ 2
m
1/ 2
f1/ 2  f ( xm , y1/ 2 )
~
f1  f ( xr , y1 )
1/ 2
l
y1  yl  dx  f1/ 2
1/ 2
Steps in algorithm







Complications: have a vector of functions and vector
of yl values
call derivs(xl,yl,f0,nord) (sets all function f0 values)
~
do n=1,nord
y1/ 2 (n)  yl (n)  dxh f 0 (n)
~
~
~
call derivs(xm, y1/ 2 , f1/ 2,nord) (sets function f1/ 2 values)
~
do n=1,nord
y1/ 2 (n)  yl (n)  dxh f1/ 2 (n)
call derivs(xm, y1/ 2 , f1/ 2,nord) (sets function f1/ 2 values)
y1 (n)  yl (n)  dx  f1/ 2 (n)
do n=1,nord
~
~
 call derivs(xr, y1 , f1 ,nord) (sets function f1 values)

Calculate R-K formula & derivs
subroutine
~
~
dx
do i=1,nord y (n)  y (n)   f (n)  2 f (n)  f (n)  f (n)
6
r
l
0
1/ 2

That’s the end of the R-K routine

subroutine derivs(x,y,yp,nord)
do n=1,nord-1:
yp(n)=y(n+1)
Lastly set yp(nord)=g(x,y,y´,…,y(nord-1))


1/ 2
1
Summary




There are a few issues with the adaptive step size
algorithm you need to be concerned about (avoiding
divide by zero etc.)
Second order systems can be turned into coupled first
order systems
Nth order systems can be turned into n coupled first
order systems
The adaptive R-K algorithm for vectors is in principle
similar to that for a single function

However, must loop over vectors
NON-EXAMINABLE BUT USEFUL TO KNOW
Implicit Methods: Backward Euler
method

Recall in the Forward Euler method
yn 1  yn
yn 1  yn  hf (t n , yn ) 
 f (tn , yn )
h

Rather than predicting forward, we can predict backward using
the value of f(tn,yn) yn  yn 1
h


 f (t n , yn )
Rewriting in terms of n+1 and n
yn 1  yn
 f (t n 1 , yn 1 )  yn 1  yn  hf (t n 1 , yn 1 )
h
Replacing yn+1=r we need to use a root finding method to solve
r  yn  hf (tn1 , r )  0
Notes on implicit Methods

Implicit methods tend to be more stable, and for the same step
size more accurate

The trade-off is the increased expense of using an interative procedure to
find the roots


Richardson Extrapolation can also be applied to implicit ODEs,
results in higher order schemes with good convergence properties


This can become very expensive in coupled systems
Use exactly the same procedure as outline in previous lecture, compare
expansions at h and h/2
Crank-Nicholson is another popular implicit method


Relies upon the derivative at both the start and end points of the interval
Second order accurate solution yn  yn 1 1
h

2
 f (tn1 , yn1 )  f (tn , yn )
Multistep methods



Thus far we’ve consider self starting methods that use
only values from xn and xn+1
Alternatively, accuracy can be improved by using a
linear combination of additional points
Utilize yn-s+1,yn-s+2,…,yn to construct approximations to
derivatives of order up to s, at tn

Example for s=2
y 'n  f (t n , yn )
f (t n , yn )  f (t n 1 , yn 1 )
y ' 'n 
h
Comparison of single and multistep
methods
Second order method

We can utilize this relationship to describe a multistep second
order method
h2
yn 1  yn  hy 'n 
 yn 

s=1 recovers the Euler method
Implicit methodologies are possible as well


h
3 f (tn , yn )  f (tn1 , yn1 )
2
Generalized to higher orders, these methods are known as
Adams-Bashforth (predictor) methods


2
y ' 'n
Adams-Moulton (predictor-corrector) methods
Since these methods rely on multisteps the first few values of y
must be calculated by another method, e.g. RK

Starting method needs to be as accurate as the multistep method
“Stiff” problems

Definition of stiffness:

“Loosely speaking, the initial value problem is referred to as being stiff if
the absolute stability requirement dictates a much smaller time step than
is needed to satisfy approximation requirements alone.”



Fomally: An IVP is stiff in some interval [0,b] if the step size needed to
maintain stability of the forward Euler method is much smaller than the step
size required to represent the solution accurately.
Stability requirements are overriding accuracy requirements
Why does this happen?

Trying to integrate smooth solutions that are surrounded by strongly
divergent or oscillatory solutions

Small deviations away from the true solution lead to forward terms being very
inaccurate
Example



Consider: y’(t)=-15y(t),
t≥0, y(0)=1
Exact solution: y(t)=e-15t,
so y(t)→0 as t→0
If we examine the
forward Euler method,
strong oscillatory
behaviour forces us to
take very small steps
even though the function
looks quite smooth
Implicit methods in stiff problems


Because implicit methods can use longer timesteps, they
are strongly favoured in integrations of stiff systems
Consider a two-stage Adams-Moulton integrator:
h
 f (tn , yn )  f (tn1 , yn1 ) 
2
 15h 
 15h 
 1 
 yn 1  1 
 yn
2 
2 


 15h 
1

2
 yn
 yn 1  
 1  15h 


2 

yn 1  yn 
Adams-Moulton solution for h=0.125
Much better behaviour and convergence
Choosing stiff solvers

This isn’t as easy as you might think



Performance of different algorithms can be quite dependent
upon the specific problem
Researchers often write papers comparing the
performance of different solvers on a given problem
and then advise on which one to use
This is a sensible way to do things


I recommend you do the same if you have a stiff problem
Try solvers from library packages like ODEPACK or the
Numerical recipes routines
Stability of the Forward Euler
method

Stability is more important than the truncation error



Consider y’=ly for some complex l
Provided Re l < 0 solution is bounded
Substitute into the Euler method
yn1  (1  lDt ) yn  (1  lDt ) 2 yn1  (1  lDt ) n1 y0

For yn to remain bounded we must have
| 1  lDt | 1

Thus a poorly chosen Dt that breaks the above inequality will
lead to yn increasing without limit
Dt 
2
|l |
Behaviour of small errors


We considered the previous y’=ly equation because it
describes the behaviour of small changes
Suppose we have a solution ys=y+e where e is the small
error.

Substitute into y’=f(t,y) and use a Taylor expansion
y 's  f (t , y s )
 e '  f (t , y  e )
 e '  f (t , y )  e
e '  e
df
dy
df
 O(e 2 )
dy
To leading order in e
Next lecture

Monte Carlo methods