7.6 Solving Multiple Nonlinear Algebraic Equations
197
where J (x i ) is just another notation for ∇F (x i ). The Eq. (7.12) is a linear system
with coefficient matrix J and right-hand side vector F (x i ). We therefore write this
system in the more familiar form
J (x i )δ = −F (x i ),
where we have introduced a symbol δ for the unknown vector x i+1 − x i that
multiplies the Jacobian J .
The i-th iteration of Newton’s method for systems of algebraic equations consists
of two steps:
1. Solve the linear system J (x i )δ = −F (x i ) with respect to δ.
2. Set x i+1 = x i + δ.
Solving systems of linear equations must make use of appropriate software. Gaussian elimination is the most common, and in general the most robust, method for this
purpose. Python’s numpy package has a module linalg that interfaces the wellknown LAPACK package with high-quality and very well tested subroutines for
linear algebra. The statement x = numpy.linalg.solve(A, b) solves a system
Ax = b with a LAPACK method based on Gaussian elimination.
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) faster methods are based on iterative techniques.
7.6.4 Implementation
Here is a very simple implementation of Newton’s method for systems of nonlinear
algebraic equations:
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
197
where J (x i ) is just another notation for ∇F (x i ). The Eq. (7.12) is a linear system
with coefficient matrix J and right-hand side vector F (x i ). We therefore write this
system in the more familiar form
J (x i )δ = −F (x i ),
where we have introduced a symbol δ for the unknown vector x i+1 − x i that
multiplies the Jacobian J .
The i-th iteration of Newton’s method for systems of algebraic equations consists
of two steps:
1. Solve the linear system J (x i )δ = −F (x i ) with respect to δ.
2. Set x i+1 = x i + δ.
Solving systems of linear equations must make use of appropriate software. Gaussian elimination is the most common, and in general the most robust, method for this
purpose. Python’s numpy package has a module linalg that interfaces the wellknown LAPACK package with high-quality and very well tested subroutines for
linear algebra. The statement x = numpy.linalg.solve(A, b) solves a system
Ax = b with a LAPACK method based on Gaussian elimination.
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) faster methods are based on iterative techniques.
7.6.4 Implementation
Here is a very simple implementation of Newton’s method for systems of nonlinear
algebraic equations:
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