Page - 197 - in Programming for Computations – Python - A Gentle Introduction to Numerical Simulations with Python 3.6, Volume Second Edition

7.6 SolvingMultipleNonlinearAlgebraicEquations 197 whereJ(xi) is just another notation for∇F(xi). The Eq.(7.12) is a linear system with coefficientmatrixJ and right-handside vectorF(xi). We thereforewrite this systemin the morefamiliar form J(xi)δ=−F(xi), where we have introduced a symbol δ for the unknown vector xi+1 − xi that multiplies theJacobianJ. The i-thiterationofNewton’smethodforsystemsofalgebraicequationsconsists of two steps: 1. Solve the linear systemJ(xi)δ=−F(xi)with respect toδ. 2. Setxi+1 =xi+δ. Solving systems of linear equations must make use of appropriate software. Gaus- sianeliminationis themostcommon,andingeneralthemostrobust,methodforthis purpose. Python’s numpy package has a module linalg that interfaces the well- known LAPACK package with high-quality and very well tested subroutines for linear algebra. The statementx = numpy.linalg.solve(A, b) solves a system Ax=bwith aLAPACK methodbasedonGaussianelimination. When nonlinear systems of algebraic equations arise from discretization of partial differential equations, the Jacobian is very often sparse, i.e., most of its elements are zero. In such cases it is important to use algorithms that can take advantage of the many zeros. Gaussian elimination is then a slow method, and (much)fastermethodsarebasedon iterative techniques. 7.6.4 Implementation Here is averysimple 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. """ F_value = F(x) F_norm = np.linalg.norm(F_value, ord=2) # l2 norm of vector iteration_counter = 0 while abs(F_norm) > eps and iteration_counter < 100: delta = np.linalg.solve(J(x), -F_value) x = x + delta F_value = F(x) F_norm = np.linalg.norm(F_value, ord=2) iteration_counter = iteration_counter + 1