736 Chapter 16.Integration of Ordinary Differential Equations Of course,we give up accuracy in following the evolution towards equilibrium if we use large stepsizes,but we maintain stability. These considerations can easily be generalized to sets of linear equations with constant coefficients: y=-C·y (16.6.9) where C is a positive definite matrix.Explicit differencing gives yn+1=(1-Ch)·yn (16.6.10) Now a matrix An tends to zero as n-oo only if the largest eigenvalue of A has magnitude less than unity.Thus y is bounded as n-oo only if the largest 菲 eigenvalue of 1-Ch is less than 1,or in other words 爱州 ICAL 2 h< (16.6.11) 入max where Amax is the largest eigenvalue of C. On the other hand,implicit differencing gives yn+1 =yn +hyn+1 (16.6.12) or 9 yn+1=(1+Ch)-1·yn (16.6.13) 9 If the eigenvalues of C are A,then the eigenvalues of(1+Ch)-1 are (1+h)-1 which has magnitude less than one for all h.(Recall that all the eigenvalues of a IENTIFIC( positive definite matrix are nonnegative.)Thus the method is stable for all stepsizes 61 h.The penalty we pay for this stability is that we are required to invert a matrix at each step. Not all equations are linear with constant coefficients,unfortunately!For the system y=f(y) (16.6.14) Numerical 10521 implicit differencing gives 43106 yn+1 yn+hf(yn+1) (16.6.15) In general this is some nasty set of nonlinear equations that has to be solved iteratively at each step.Suppose we try linearizing the equations,as in Newton's method: Yn+1=yn +h + (16.6.16) Here of/Oy is the matrix of the partial derivatives of the right-hand side(the Jacobian matrix).Rearrange equation(16.6.16)into the form Yn+1=yn+h 1-hof] ·fyn) (16.6.17736 Chapter 16. Integration of Ordinary Differential Equations Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machine￾readable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). Of course, we give up accuracy in following the evolution towards equilibrium if we use large stepsizes, but we maintain stability. These considerations can easily be generalized to sets of linear equations with constant coefficients: y
= −C · y (16.6.9) where C is a positive definite matrix. Explicit differencing gives yn+1 = (1 − Ch) · yn (16.6.10) Now a matrix An tends to zero as n → ∞ only if the largest eigenvalue of A has magnitude less than unity. Thus yn is bounded as n → ∞ only if the largest eigenvalue of 1 − Ch is less than 1, or in other words h < 2 λmax (16.6.11) where λmax is the largest eigenvalue of C. On the other hand, implicit differencing gives yn+1 = yn + hy
n+1 (16.6.12) or yn+1 = (1 + Ch) −1 · yn (16.6.13) If the eigenvalues of C are λ, then the eigenvalues of (1 + Ch) −1 are (1 + λh)−1, which has magnitude less than one for all h. (Recall that all the eigenvalues of a positive definite matrix are nonnegative.) Thus the method is stable for all stepsizes h. The penalty we pay for this stability is that we are required to invert a matrix at each step. Not all equations are linear with constant coefficients, unfortunately! For the system y
= f(y) (16.6.14) implicit differencing gives yn+1 = yn + hf(yn+1) (16.6.15) In general this is some nasty set of nonlinear equations that has to be solved iteratively at each step. Suppose we try linearizing the equations, as in Newton’s method: yn+1 = yn + h
f(yn) + ∂f ∂y yn · (yn+1 − yn) (16.6.16) Here ∂f/∂y is the matrix of the partial derivatives of the right-hand side (the Jacobian matrix). Rearrange equation (16.6.16) into the form yn+1 = yn + h  1 − h ∂f ∂y −1 · f(yn) (16.6.17)