APC523AST523 Scientific Computation in Astrophysics

1
APC523/AST523Scientific Computation in
Astrophysics

• Lecture 7 (cont.)
• Nonlinear Equations

2
Nonlinear systems
Perhaps the next easiest problem to solve in
numerical analysis (after linear algebra) is the
solution of a nonlinear equation in a single
variable. Amounts to finding the roots of
f(x) 0 (Any nonlinear equation can be cast in
this form.) Solving a system of nonlinear
equations is MUCH harder. Hardly anything that
works for one equation in one unknown works for
multi-variable systems. In fact, cannot even
guarantee a solution exists in the multi-variable
problem. Discuss later.
3
What sort of problems lead to nonlinear equations?

• A variety of physical effects are modeled with
  nonlinear equations, for example
• Calculation of thermal equilibrium in plasma
  (balance of heating and cooling rates, the latter
  is strongly nonlinear function of T)
• Equilibrium reaction rates in interstellar
  chemistry (quadratic in abundances of each
  species)
• Nuclear reaction networks relevant in stellar
  interiors, explosions, etc. (quadratic in number
  density of species, highly nonlinear function of
  T)
• In addition, some basic numerical algorithms lead
  to nonlinear equations, e.g.
• Finding eigenvalues of a matrix is equivalent to
  finding the roots of f(l) det(lI - A)
• Iterative methods for solving sparse systems
  (e.g. conjugate gradient) involve root finding
• Iterative methods for function evaluation are
  basically root finding (e.g. evaluating sqrt(x))

4
Constraints on methods
If evaluating function f is expensive, want to
find method that minimizes function
evaluations. If the derivative of f is known,
this can be used to build faster methods.
5
It can be very helpful (important) to know
something about what the function f(x) looks
like. Plot it.
ophirgt sm Hello James, please give me a command
set x1,100 set x0.1x - 5. set
y1.-2./(xx1.) device x11 limits x y
box connect x y
6
How to find roots.
First, it is helpful to know that the root is
bracketed. Consider an interval a,b. If
f(a)f(b) lt 0, there must be at least one root in
the interval, if f(x) is continuous.
7
But how to find the interval a,b? No good
method in general. Simplest method is often the
best. Start at point a. Compute f(a). Take
increasingly larger steps away until f(b) changes
sign. But beware!
8
Roots may be so closely spaced it is easy to miss
them.
Or there may be more than one (but an odd number)
of roots in the domain.
9
Bisection.
Once the root has been bracketed, the simplest
way to converge on it is bisection. Compute
midpoint c (ab)/2. Evaluate f(c) If f(c)0,
then root has been found, and we are
done. Otherwise, choose new domain a,c if
f(a)f(c) lt 0
c,b if f(c)f(b) lt 0
10
Graphical demonstration of bisection.
a c b
11
Convergence of bisection.
Bisection always works. Converges linearly. To
reach a given error tolerence e, it takes
log2(e/e0) steps. However, if domain contains
more than one root, which root the method
converges on depends on the initial choice of
a,b Practical error tolerence is important. If
round-off error in evaluating f(x) suffers
amplification by k, then will never be able to
converge to machine precision.
12
Secant method
Can we do better than linear convergence? Smooth
functions are always locally linear approximate
with a straight line. That is, draw a straight
line between a and b, root will lie near the
intersection of this line with the x-axis.
The formula is
13
Convergence of secant method.
Convergence rate is faster than linear
(superlinear), but not uniform. For example, for
f(x) x2 - 2 in a,b 1,2 Converges to
10 digits of accuracy in 6 iterations. This is
the same iterative method for finding sqrt
introduced in lecture 5.
14
But secant can also diverge.
3
2
1
Would not happen if root is kept bracketed, that
is use 3 and 1 instead of 3 and 2. This is the
false position method. Converges slower, but
is more sure.
15
Brents method
Probably the best method available if only the
function, but not its derivative, is
available. Combines bisection, root bracketing,
and quadratic convergence. Guaranteed to
converge on root provided function can be
evaluated anywhere in interval containing
root. Fits a quadratic function to the three
previous points at each iteration. Next guess is
taken to be zero crossing of this function. If
interpolation function will result in a new guess
that lies outside interval, a bisection step is
taken instead. See NR, section 9.3
16
Newton-Raphson iteration
Faster convergence is possible if the derivative
of f is available. Expand f(x) in a Taylor
series d is the correction added to
current guess of root i.e. xi1 xi
d(xi)
17
Graphically, NR uses tangent line at current
guess to correct guess. Zero crossing of tangent
line is next guess.
18
Failure of NR
NR only works in neighborhood of root d/x ltlt 1.
If slope is close to zero, NR fails
spectacularly. If second derivative changes sign
at root, NR may never converge.
But generally NR extremely robust, especially
when combined with bracketing. Get quadratic
convergence.
19
Astrophysical exampleCooling in the
interstellar medium
20
Calculating the temperature

• Recall the ISM gas and dust in the space
  between stars.
• Question What is the equilibrium temperature of
  the ISM?
• Must balance the net heating and cooling rates.
• Heating is produced by
• photons from stars
• cosmic rays (high energy particles)
• Both processes are essentially independent of the
  T of the ISM
• Cooling is produced by
• Bremsstrahlung photons produced by collisions
  between ions and electrons
• Inelastic collisions followed by radiative decay.
  Many, many species contribute to the net
  process, each important at a different T.
• Thermal radiation by dust heated by collisions
  with gas
• These processes are all nonlinear functions of T.

21
Per-particle cooling rate is a nonlinear function
of T. Balance of H C(T) can be found
graphically. However, if many solutions
needed for many different H, may need
computational methods.
H
22
APC523/AST523Scientific Computation in
Astrophysics

• Lecture 8
• Nonlinear Equations (cont.)

23
Systems of nonlinear equations
Solving a system of nonlinear equations is MUCH
harder. Hardly anything that works well for
one equation in one unknown works well for
multi-variable systems of equations. In fact,
cannot even guarantee a solution exists in the
multi-variable problem. Consider F(x,y)
0 G(x,y) 0 Solutions are the x,y values that
satisfy both equations. But there may be NONE.
Moreover, difficult to tell without solving both
equations simultaneously and looking for
overlap. Most methods are generalizations of 1D
methods.
24
Multidimensional NR
If you are already near a root, NR is perhaps the
best choice, since it converges
quadratically. In one variable, NR can be
thought of as a method to find where the function
g(x) x, where The multidimensional extension
of this is solving G(x) x, where J is the
Jacobian of F(x). Of course, computing the
inverse J-1 is prohibitive,
25

• Instead, use two-step NR iteration
• Starting with an initial guess for the solution
  x0
• Solve J(x0)y -F(x0) for y
• Set x1 x1 y
• If y lt e STOP. Else repeat using x1.
• Note Step (1) requires solving a NxN matrix
  requires methods for linear algebra (LU
  decomposition).

Converges quadratically, IF the initial guess for
the solution is close the the root. Very costly.
Requires computing NxN matrix every iteration,
and O(N3) operations to solve step (1). Matrix J
may be very difficult to compute (singular, or
analytic formulae for derivatives not known).
26
Multidimensional Secant Method (Broydens Method)
Cheaper than NR, but does not converge as fast
except very near to root. Basic idea is to (1)
replace evaluation of Jacobian with
finite-difference approximation that only
requires F(x) to be evaluated, and (2) find an
approximation for J-1 at each iteration, using
value at last iteration (Sherman-Morrison
Formula). Both are O(N2) operations. Then
solving J(x0)y -F(x0) for y requires only
matrix-vector product y -J-1(x0
)F(x0). This is O(N2) operations. So, Broydens
method can be much faster. Reduces to Secant
Method in one-dimension Note both NR and
Broydens method fail when J is singular!
27
Globally convergent methods Steepest descent
If a good initial guess for the root is NOT
available, then need a globally convergent
method, e.g. steepest descent is useful (but
slow). Also very useful as a way to get a good
first guess for NR or Broydens method. The idea
is to find the minimum of the single
function Newtons method follows the descent
direction of G However, taking a full Newton
step can overshoot root when quadratic
approximation fails (far from root). So we want
to limit Newton step by some parameter a, where 0
lt a lt 1. Various techniques for finding a, based
on magnitude of grad(G) at xold and xnew.
28
Example Interstellar Chemistry Luminous blue
stars which must be young are always found
near interstellar clouds Conclusion stars
form in dense interstellar gas clouds.
29
Eagle nebula (M16)
30
Images of Eagle nebula at radio wavelengths
reveal kinematics of molecular gas
31
BU-FCRAO Galactic Ring Survey
13CO column density
32
Why is chemistry important?
Most of the gas is molecular. Dominant species
is H2, but complex molecules like PAHs
important. Only trace species (CO, NH3 , ) are
detectable. Interpreting spectra requires
understanding abundances. Molecular emission
line radiation controls the temperature. Chemical
reactions control the ionization fraction (and
therefore coupling of gas to interstellar
magnetic fields).
33
PAH emission features from a reflection nebulae
Structure of 4 PAHs
34
How to use chemistry to compute temperature
needed for dynamics.

• Need to know molecular abundances to compute the
  total (per volume) cooling rate L(T) due to line
  emission.
• Use equilibrium chemical model to compute
  abundances (assuming reaction times short
  compared to dynamical times -- may not be true in
  shocks).
• Requires knowing photoionization rate, cosmic ray
  ionization rate, and then solving chemical model
  for 100s of species.
• Reaction rates poorly known
• Ionization rates uncertain
• Temperature comes from balancing heating versus
  cooling (including mechanical heating and/or
  cooling terms).

35
Current state-of-the-art.

• Time-dependent chemical models embedded in 3D MHD
  dynamics is prohibitive.
• Instead Time-dependent chemistry with no
  dynamics
• Even simple astrochemical problems remain to be
  understood
• How does H2 form on grains?
• Abundance of H3
• Diffuse interstellar bands (unidentified
  absorption lines).

36
Diffuse interstellar bands -- unidentified
spectral features Ultra-large molecules or
ultra-small grains?
37
Elements of equilibrium chemical reaction
networks Consider reaction between two species
A and B A B -gt AB rate nA nB RAB The
reverse is also possible AB -gt A B, rate
nAB RAB In equilibrium, these rates are equal.
nA nB RAB nAB RAB Combine with
normalizations nA nAB const. nB nAB
const. Get quadratic equation in nAB which is
easily solved. But if there are many possible
reactions involving A, B, C, . And many products
AB, AC, BC, , then problem is multidimensional
nonlinear root finding for abundances of
products. Example is UMIST Database for
astrochemistry http//www.rate99.co.uk
38
Example nuclear reaction networks in stars
39
Proton-proton chainMost important reaction in
Sun is PPI chain
1H 1H -gt 2H e n 2H 1H -gt 3He
g 3He 3He -gt 4He 1H 1H Net result is 4H
fused into 4He
e positron n neutrino g photon
40

• 69 of the time, H fusion occurs via PPI chain in
  Sun
• 31 of the time, PPII chain occurs
• 1H 1H -gt 2H e n
• 2H 1H -gt 3He g
• 3He 4He -gt 7Be g
• 7Be e- -gt 7Li n
• 7Li 1H -gt 2 4He
• Other chains (CNO cycle) are negligible, but
  important in other stars

41
PPI Chainstars with M lt 2MSun
42
CNO cycle stars with M gt 2MSun
Uses C, N, and O nuclei to catalyze fusion of H
into He
43
Elements of equilibrium nuclear reaction
networks Consider reaction between two species
i and j i j -gt k l rate ljk In
general, the time rate of change of the mole
fraction of species i obeys a differential
equation. Single body reactions (beta-decay),
and three-body reactions (triple-a reaction) can
also be included as extra terms. Ignore
diffusion of Yi, and assume equilibrium (dY/dt
0). Then RHS is system of nonlinear equations
for the mole fractions Yi. However, new
normalizations are needed compared to chemistry.
Number of nuclei of i,j not conserved. Instead,
conserve lepton number, energy, etc. Often, the
system of rate equations must be supplemented
with other equations for the temperature and
density (assuming pressure fixed). In the case
of stellar structure, ODEs must be solved to
compute r, P, T.
44
Detection of solar neutrinos produced by p-p chain

• Underground tank of C2Cl4
• Detects 37Ar produced by reaction between 37Cl
  and n
• Sensitive to production of few Ar atoms per month

45
Solar neutrino results

• For 30 years, measured rate was 1/3 predicted
  rate from standard solar model solar neutrino
  problem
• In last 5 years, it has been discovered that
  neutrinos have mass undergo flavor mixing
• New experiments (e.g. SNO) sensitive to all
  flavors of neutrinos find rates consistent with
  standard solar model
• The problem was not with the solar model, but
  our understanding of particle physics!

46
APC523/AST523Scientific Computation in
Astrophysics

• Lecture 8
• ODEs I Initial Value Problems

47
Ordinary differential equations
Involve derivative wrt one independent
variable. As well see, huge number of problems
in astrophysics result in ODEs. Problems
involving ODEs can always be reduced to a system
of first-order equations. For example Can be
rewritten as so only need to consider
methods for systems of first order equations
48
Boundary conditions

• Must also specify boundary conditions to find a
  unique solution.
• Boundary conditions are algebraic conditions on
  the values of yi at discrete point(s) xi.
  Systems of first order equations require only one
  BC per equation.
• Two categories of BCs
• Initial value problems all yi specified at same
  starting point x0, and the solution is integrated
  from x0 to final point xf
• Boundary value problems yi are specified at two
  or more points, e.g. some are specified at x0 and
  some at xf (two-point BV problem).
• Latter is extremely important in astrophysics
  equations of stellar structure
• IVPs are generally easier to solve, so start with
  them.

49
Eulers method
Consider the first order ODE Approximate the
derivative with a forward finite-difference Then
, starting from a known value yn at xn (the
initial value), the solution is Here h is the
stepsize, and the xn are the mesh
points. Integration is unsymmetric. Relies only
on information at xn. Cannot be reversed (this
is a big problem!). Method is first-order
accurate (leading error is O(h)).
50
Runge-Kutta methods
Recall we symmetrized derivatives (centered
difference formula) to get higher-order
accuracy. Can do the same here, that is compute
the derivative f at the midpoint of the step,
rather than at xn. To get values at midpoint,
use first-order Euler. So the method
becomes Can show leading error is O(h2), so
this is 2nd-order RK.
51
Fourth-order RK
Recall we used Richardson extrapolation to
combine quadratures with different stepsizes with
different weights so as to cancel the leading
error term. Conceptually, can do the same thing
here. Higher-order RK methods result from
combining lower-order methods so as to cancel
leading error term. Classic is 4th-order RK.
Let Then Biggest disadvantage of method is
it requires f to be evaluated 4 times per
integration step.
52
RK with adaptive stepsizes
Ideally we want to vary stepsize h according to
magnitude of f. When f big, take small
steps When f small, take bigger steps
y
x
small steps
53
Cost saving with adaptive steps
Can take orders of magnitude fewer steps with
adaptive stepping. But how much overhead is
added to the method? Again, basic idea is to use
something like Richardson extrapolation. Take
one step at 2h (costs 4 evaluations of f) Take
two steps at h (costs 8-17 evaluations of
f) Without adaptive steps, 8 evaluations of f
would be needed (two steps at h). Remember,
were trying to make the steps bigger. So
overhead is (47)/8 1.375 Actually very little
overhead is added!
54
Step doubling
For fourth-order RK, the leading term in the
truncation error is O(h5) Suppose take steps of
size h1 ---gt will get error D1 Suppose take steps
of size h0 ---gt will get error D0 Then
Rewrite as
(1) Now let D0 be desired accuracy of
solution. Take a step of size h1. Estimate D1.
Then (1) gives stepsize necessary to achieve
desired accuracy.
55
But making stepsize smaller leads to larger
accumulated error for fixed D0. So we want to
make D0 smaller as h gets smaller to keep
accumulated error constant. If we assume D
proportional to h, we get
(2) How to choose
between (1) and (2). Take pragmatic
approach When D0 gt D1 use D0 lt D1
use This gives smallest rate of increase of
stepsize h largest rate of
decrease of stepsize h
56
Embedded RK
Another way to estimate the local truncation
error for adaptive stepsize methods. RK methods
of order M require M evaluations of f, for M lt
5. Thus the popularity of RK4 --- highest order
with only M evaluations of f Fehlberg found an
RK5 formula that requires only 6 evaluations of
f. Moreover, these same 6 values could be used
to construct a RK4 algorithm. Comparing the RK4
and RK5 methods allows an estimate of the local
truncation error in the RK4 integrator, and a way
to control adaptive stepsize ---gt embedded RK
methods.
57
Modified midpoint method
RK and Euler methods require data from mesh point
n to construct solution at n1. Wouldnt it be
better to use a centered difference form to
estimate the derivative Problems are
(1) how to start the method? A Use forward
Euler for first step. (2) how to end the
method? A Use modified midpoint to get yN
from yN-1 and yN-2. Use forward Euler to get yN
from yN-1. Then average two. Why use forward
Euler? Makes error at start and end of domain
symmetric. Get 2nd order method with only one
evaluation of f. Cheaper than RK2. Its
reversible! Related to leap-frog scheme.
58
But the real power of modified midpoint is the
result due to Gragg. The truncation error over
the interval x,xH contains only even powers of
h. So by using Richardson extrapolation, can
increase order of convergence by two powers of h
each time. Thus is fourth-order
accurate. Just as in Romberg integration
methods, can extend this idea to create higher
order methods, leading to
59
Bulirsch-Stoer method

• Probably most efficient method for high accuracy
  integration of smooth functions.
• Based on three ideas. AT EACH INTEGRATION STEP
• Use modified midpoint to integrate equations, so
  that truncation error is an even power of h.
  Perform multiple integrations using decreasing
  stepsize that follows sequence hH/n
  n2,4,6,8,10,12,14,16,
• Use Richardson extrapolation to measure how
  truncation error converges with decreasing
  stepsize h.
• Fit rational functions (R P/Q these behave
  better at poles than polynomials) to fit
  convergence of error with h. Use this functional
  fit to estimate solution for h0 (analytic
  solution).
• In practice, polynomial fitting is often used
  instead of rational functions.
• Functional fitting in step (3) also provides
  error estimate. Can repeat integration of step
  (1) with ever higher orders until error satisfies
  criterion. Or, can adapt stepsize H.

60
Multistep methods
Euler, RK, and Bulirsch-Stoer are all single step
methods, that is they only require yn at xn to
integrate to next step. (Although inside BS is
the modified midpoint method, which is
multistep.) Multistep methods fit polynomial of
degree M-1 to yn-i i0,..,M data
points. Probably best known is fourth-order
Adams-Bashforth Has to be started with
single-step method. Difficult to implement
adaptive stepsizes, since interpolation function
assumes constant h.
61
Predictor-corrector
A large class of methods designed to get
high-order accuracy for smooth functions.
Probably the most popular example is
Adams-Bashforth-Moulton. Adams-Bashforth is
one-sided only uses data from xn and previous
values to estimate xn1 Better is to use
symmetric method to estimate f(xn, yn). Use
3rd-order Adams-Bashforth to estimate xn1. This
is predict step. Then fit polynomial to n-1, n,
n1 to estimate f(xn, yn). Use to recompute
yn1. This is correct step. Large number of
legacy codes use predictor-corrector methods.
But Bulirsch-Stoer is probably easier to code and
just as accurate.