Nonlinear EquationsNonlinear Equations●Often (most of the time??) the relevant system of equations is not linear in the unknowns.●Then, cannot decompose as Ax = b. Oh well.●Instead write as:(1) f(x) = 0 function of one variable (1-D)(2) f(x) = 0 x = (x1,x2,...,xn), f = (f1,f2,...,fn) (n-D)●Not guaranteed to have any real solutions, but generally do for astrophysical problems.Solutions in 1-DSolutions in 1-D●Generally, solving multi-D equations is much harder, so we'll start with the 1-D case first...●By writing f(x) = 0 we have reduced the problem to solving for the roots of f.●In 1-D it is usually possible to trap or bracket the desired roots and hunt them down.●Typically all root-finding methods proceed by iteration, improving a trial solution until some convergence criterion is satisfied.Function PathologiesFunction Pathologies●Before blindly applying a root-finding algorithm to a problem, it is essential that the form of the equation in question be understood: graph it!●For smooth functions good algorithms will always converge, provided the initial guess is good enough.●Pathologies include discontinuities, singularities, multiple or very close roots, or no roots at all!Numerical Root FindingNumerical Root Finding●Suppose f(a) and f(b) have opposite sign.●Then, if f is continuous on the interval (a,b), there must be at least one root between a and b (this is the Intermediate Value Theorem).●Such roots are said to be bracketed.Bracketed rootMany rootsNo rootsExample ApplicationExample Application●Use root finding to calculate the equilibrium temperature of the ISM.●The ISM is a very diffuse plasma.–Heated by nearby stars and cosmic rays.–Cooled by a variety of processes:●Bremsstrahlung: collisions between electrons and ions●Atom-electron collisions followed by radiative decay●Thermal radiation from dust grainsExample, Cont'dExample, Cont'd●Equilibrium temperature given when: Rate of Heating H = Rate of Cooling C●Often H is not a function of temperature T.●Usually C is a complex, nonlinear function of T.●To solve, find T such that H - C(T) = 0.Bracketing and BisectionBracketing and Bisection●NRiC 9.1 lists some simple bracketing routines.●Once bracketed, root is easy to find by bisection:–Evaluate f at interval midpoint (a + b) / 2.–Root must be bracketed by midpoint and whichever a or b gives f of opposite sign.–Bracketing interval decreases by 2 each iteration: εn+1 = εn / 2.–Hence to achieve error tolerance of ε starting from interval of size ε0 requires n = log2(ε0/ε) steps.ConvergenceConvergence●Bisection converges linearly (first power of ε).●Methods in which εn+1 = (constant) × (εn)mm > 1 are said to converge superlinearly.●Note error actually decreases exponentially for bisection. It converges "linearly" because successive figures are won linearly with computational effort.ToleranceTolerance●What is a practical tolerance ε for convergence?●Cannot be less than round-off error.●For single-precision (float) accuracy, typically take ε = 10-6 in fractional error. i.e.where x = numerical solution, xr = actual root.●When f(xr) = 0 this fails, so use ε = 10-6 as absolute error (or perhaps use ε(|a| + |b|)/2).  −     ~ −Newton-Raphson MethodNewton-Raphson Method●Can one do better than linear convergence? Duh!●Expand f(x) in a Taylor series: f(x + δ) = f(x) + f '(x) δ + f "(x) δ2/2 + ...●For δ << x, drop higher order terms, so: f(x + δ) = 0 ⇒ δ = - f(x) / f '(x)●δ is correction added to current guess of root: i.e. xi+1 = xi + δNewton-Raphson, Cont'dNewton-Raphson, Cont'd●Graphically, Newton-Raphson (NR) uses tangent line at xi to find zero crossing, then uses x at zero crossing as next guess:●Note: only works near root (δ << x)...Newton-Raphson, Cont'dNewton-Raphson, Cont'd●When higher order terms important, NR fails spectacularly. Other pathologies exist too:Shoots to infinityNever convergesNewton-Raphson, Cont'dNewton-Raphson, Cont'd●Why use NR if it fails so badly?●Rate of convergence: εi+1 = εi - f(xi) / f '(xi)●Taylor expand f(xi) & f '(xi) to get: εi+1 = - εi2 f ''(xi) / f '(xi) [quadratic!]●Note both f(x) and f '(x) must be evaluated each iteration, plus both must be continuous near root.●Best use of NR is to "polish-up" bisection root.Nonlinear Systems of EquationsNonlinear Systems of Equations●Consider the system f(x,y) = 0, g(x,y) = 0. Plot zero contours of f & g:●No information in f about g, and vice versa.–In general, no good method for finding roots.Nonlinear Systems, Cont'dNonlinear Systems, Cont'd●If you are near root, best bet is NR.e.g. For F(x) = 0, choose xi+1 = xi + δ, where F'(x) δ = - F(x)●This is a matrix equation: F'(x) is a matrix with elements ∂Fi / ∂xj. The matrix is called the Jacobian.●Written out:∂ ∂ ∂ ∂ = −     ∂ ∂ ∂ ∂ = −    Nonlinear Systems, Cont'dNonlinear Systems, Cont'd●Given initial guess, must evaluate matrix elements and RHS, solve system for δ, and compute next iteration xi+1. Then repeat (must solve 2 × 2 linear system each time).●Essentially the non-linear system has been linearized to make it easier to work with.●NriC 9.7 discusses a global convergence strategy that combines multi-D NR with "backtracking" to improve chances of finding solutions.Example: Interstellar ChemistryExample: Interstellar Chemistry●ISM is multiphase plasma consisting of electrons, ions, atoms, and molecules.●Originally, the ISM was thought to be too hostile for molecules.●But in 1968-69, radio observations discovered absorption/emission lines of NH3, H2CO, H2O, ... ●Lots of organic molecules, e.g. CH3CH2OH (ethanol).Example, Cont'dExample, Cont'd●In some places, all atoms have been incorporated into molecules.●For example, molecular clouds: dense, cold clouds of gas composed primarily of molecules.(T ~30 K, n ~106 cm-3, M ~105-6 M⊙, R ~ 10 -100 pc).●How do you predict what abudances of different molecules should be, given n and T?●Need to solve a chemical reaction network.Example, Cont'dExample, Cont'd●Consider reaction between two species A and B: A + B → AB reaction rate = nAnBRAB●Reverse also possible: AB → A + B reaction rate = nABR'AB●In equilibrium:(1) nAnBRAB = nABR'AB(2) nA + nAB =