Page - 205 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python

6.6 SolvingMultipleNonlinearAlgebraicEquations 205 6.6.3 Newton’sMethod TheideaofNewton’smethodis thatwehavesomeapproximationxi to therootand seek a new (andhopefully better) approximationxiC1 byapproximatingF.xiC1/ by a linear function and solve the corresponding linear systemof algebraic equa- tions.Weapproximate thenonlinearproblemF.xiC1/D0by the linearproblem F.xi/CJ.xi/.xiC1 xi/D0; (6.13) whereJ.xi/ is just another notation forrF.xi/. The equation (6.13) is a linear systemwith coefficientmatrixJ and right-hand side vectorF.xi/. We therefore write this systemin themore familiar form J.xi/ı D F.xi/; wherewehave introducea symbolı for theunknownvectorxiC1 xi thatmulti- plies the JacobianJ. The i-th iteration ofNewton’smethod for systems of algebraic equations con- sists of twosteps: 1. Solve the linear systemJ.xi/ı D F.xi/with respect toı. 2. SetxiC1 Dxi Cı. Solving systemsof linear equationsmustmakeuseof appropriate software. Gaus- sian elimination is themost common, and in general themost robust, method for this purpose. Python’s numpy package has a module linalg that interfaces the well-knownLAPACKpackagewith high-quality andverywell tested subroutines for linearalgebra. Thestatementx = numpy.linalg.solve(A, b)solvesasys- temAxDbwithaLAPACKmethodbasedonGaussian elimination. Whennonlinear systemsof algebraicequationsarise fromdiscretizationofpar- tialdifferentialequations, theJacobianisveryoftensparse, i.e.,mostofitselements are zero. In such cases it is important to use algorithms that can take advantageof themany zeros. Gaussian elimination is then a slowmethod, and (much) faster methodsarebasedon iterative techniques. 6.6.4 Implementation Here isaverysimple implementationofNewton’smethodforsystemsofnonlinear algebraicequations: import numpy as np def Newton_system(F, J, x, eps): """ Solve nonlinear system F=0 by Newton’s method. J is the Jacobian of F. Both F and J must be functions of x. At input, x holds the start value. The iteration continues until ||F|| < eps. """