Newton for nonlinear systems Tobin A. Driscoll Richard J. Braun
The rootfinding problem becomes much more difficult when multiple variables and equations are involved.
Particular problems are often posed using scalar variables and equations.
The steady state of interactions between the population w ( t ) w(t) w ( t ) h ( t ) h(t) h ( t )
a h − b h w = 0 , − c w + d w h = 0 , \begin{split}
ah - b h w &= 0, \\
-cw + d w h &= 0,
\end{split} ah − bh w − c w + d w h = 0 , = 0 , for positive parameters a , b , c , d a,b,c,d a , b , c , d (4.5.1) x = [ h , w ] \mathbf{x}=[h,w] x = [ h , w ] f 1 ( x 1 , x 2 ) = a x 1 − b x 1 x 2 f_1(x_1,x_2) = ax_1 - bx_1x_2 f 1 ( x 1 , x 2 ) = a x 1 − b x 1 x 2 f 2 ( x 1 , x 2 ) = − c x 2 + d x 1 x 2 f_2(x_1,x_2)= -c x_2 + d x_1 x_2 f 2 ( x 1 , x 2 ) = − c x 2 + d x 1 x 2
While the equations of Example 4.5.1
4.5.1 Linear model ¶ To extend rootfinding methods to systems, we will keep to the basic philosophy of constructing easily managed models of the exact function. As usual, the starting point is a linear model. We first need to define what it means to take a derivative of a vector-valued function of a vector variable.
A multidimensional Taylor series begins with the linear approximation
f ( x + h ) = f ( x ) + J ( x ) h + O ( ∥ h ∥ 2 ) , \mathbf{f}(\mathbf{x}+\mathbf{h}) = \mathbf{f}(\mathbf{x}) + \mathbf{J}(\mathbf{x})\mathbf{h} + O(\| \mathbf{h} \|^2), f ( x + h ) = f ( x ) + J ( x ) h + O ( ∥ h ∥ 2 ) , where J \mathbf{J} J f \mathbf{f} f x \mathbf{x} x h \mathbf{h} h x \mathbf{x} x
The terms f ( x ) + J ( x ) h \mathbf{f}(\mathbf{x})+\mathbf{J}(\mathbf{x})\mathbf{h} f ( x ) + J ( x ) h (4.5.7) f \mathbf{f} f x \mathbf{x} x O ( ∥ h ∥ 2 ) O(\| \mathbf{h} \|^2) O ( ∥ h ∥ 2 ) f \mathbf{f} f f ( x ) = A x − b \mathbf{f}(\mathbf{x})=\mathbf{A}\mathbf{x}-\mathbf{b} f ( x ) = Ax − b A \mathbf{A} A (4.5.7)
Because of the Jacobian’s role in (4.5.7) J ( x ) \mathbf{J}(\mathbf{x}) J ( x ) f ′ ( x ) \mathbf{f}{\,}'(\mathbf{x}) f ′ ( x ) x \mathbf{x} x
4.5.2 The multidimensional Newton iteration ¶ With a method in hand for constructing a linear model for the vector system f ( x ) \mathbf{f}(\mathbf{x}) f ( x ) x k \mathbf{x}_k x k h = x − x k \mathbf{h} = \mathbf{x}-\mathbf{x}_k h = x − x k (4.5.7)
f ( x ) ≈ q ( x ) = f ( x k ) + J ( x k ) ( x − x k ) . \mathbf{f}(\mathbf{x}) \approx \mathbf{q}(\mathbf{x}) = \mathbf{f}(\mathbf{x}_k) + \mathbf{J}(\mathbf{x}_k)(\mathbf{x}-\mathbf{x}_k). f ( x ) ≈ q ( x ) = f ( x k ) + J ( x k ) ( x − x k ) . We define the next iteration value x k + 1 \mathbf{x}_{k+1} x k + 1 q ( x k + 1 ) = 0 \mathbf{q}(\mathbf{x}_{k+1})=\boldsymbol{0} q ( x k + 1 ) = 0
0 = f ( x k ) + J ( x k ) ( x k + 1 − x k ) , \begin{split}
\boldsymbol{0} &= \mathbf{f}(\mathbf{x}_k) + \mathbf{J}(\mathbf{x}_k)(\mathbf{x}_{k+1}-\mathbf{x}_k),\\
\end{split} 0 = f ( x k ) + J ( x k ) ( x k + 1 − x k ) , which can be rearranged into
x k + 1 = x k − [ J ( x k ) ] − 1 f ( x k ) . \mathbf{x}_{k+1} = \mathbf{x}_k - \bigl[\mathbf{J}(\mathbf{x}_k)\bigr]^{-1} \mathbf{f}(\mathbf{x}_k). x k + 1 = x k − [ J ( x k ) ] − 1 f ( x k ) . Note that J − 1 f \mathbf{J}^{-1}\mathbf{f} J − 1 f f / f ′ f/f' f / f ′
Algorithm 4.5.1 (Multidimensional Newton’s method)
Given f \mathbf{f} f x 1 \mathbf{x}_1 x 1 k = 1 , 2 , 3 , … k=1,2,3,\ldots k = 1 , 2 , 3 , …
Compute y k = f ( x k ) \mathbf{y}_k = \mathbf{f}(\mathbf{x}_k) y k = f ( x k ) A k = f ′ ( x k ) \mathbf{A}_k=\mathbf{f\,}'(\mathbf{x}_k) A k = f ′ ( x k ) Solve the linear system A k s k = − y k \mathbf{A}_k\mathbf{s}_k = -\mathbf{y}_k A k s k = − y k Newton step s k \mathbf{s}_k s k Let x k + 1 = x k + s k \mathbf{x}_{k+1} = \mathbf{x}_k + \mathbf{s}_k x k + 1 = x k + s k An extension of our series analysis of the scalar Newton’s method shows that the vector version is also quadratically convergent in any vector norm, under suitable circumstances and when the iteration converges at all.
4.5.3 Implementation ¶ An implementation of Newton’s method for systems is given in Function 4.5.2 norm, Function 4.5.2 Function 4.3.1
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
"""
newtonsys(f, jac, x₁ [;maxiter, ftol, xtol])
Use Newton's method to find a root of a system of equations,
starting from `x₁`. The functions `f` and `jac` should return the
residual vector and the Jacobian matrix, respectively. Returns the
history of root estimates as a vector of vectors.
The optional keyword parameters set the maximum number of iterations
and the stopping tolerance for values of `f` and changes in `x`.
"""
function newtonsys(f, jac, x₁; maxiter = 40, ftol = 1e-13, xtol = 1e-13)
x = [float(x₁)]
y, J = f(x₁), jac(x₁)
Δx = Inf # for initial pass below
k = 1
while (norm(Δx) > xtol) && (norm(y) > ftol)
Δx = -(J \ y) # Newton step
push!(x, x[k] + Δx) # append to history
k += 1
y, J = f(x[k]), jac(x[k])
if k == maxiter
@warn "Maximum number of iterations reached."
break
end
end
return x
endAbout the code
The output of Function 4.5.2 float before the iteration starts.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
function x = newtonsys(f,x1)
% NEWTONSYS Newton's method for a system of equations.
% Input:
% f function that computes residual and Jacobian matrix
% x1 initial root approximation (n-vector)
% Output
% x array of approximations (one per column, last is best)
% Stopping parameters.
funtol = 1000 * eps; % stop for small || f(x) ||
xtol = 1000 * eps; % stop for small || x change ||
maxiter = 40; % stop after this many iterations
x = x1(:);
[y, J] = f(x1);
dx = Inf;
k = 1;
while (norm(dx) > xtol) && (norm(y) > funtol) && (k < maxiter)
dx = -(J \ y); % Newton step
x(:, k+1) = x(:, k) + dx;
k = k+1;
[y, J] = f(x(:, k));
end
if k==maxiter
warning('Maximum number of iterations reached.')
end
end1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
def newtonsys(f, jac, x1):
"""
newtonsys(f, jac, x1)
Use Newton's method to find a root of a system of equations, starting from x1. The
function f should return the residual vector, and the function jac should return
the Jacobian matrix. Returns root estimates as a matrix, one estimate per column.
"""
# Operating parameters.
funtol = 1000 * np.finfo(float).eps
xtol = 1000 * np.finfo(float).eps
maxiter = 40
x = np.zeros((maxiter, len(x1)))
x[0] = x1
y, J = f(x1), jac(x1)
dx = 10.0 # for initial pass below
k = 0
while (norm(dx) > xtol) and (norm(y) > funtol) and (k < maxiter):
dx = -lstsq(J, y)[0] # Newton step
x[k+1] = x[k] + dx
k = k + 1
y, J = f(x[k]), jac(x[k])
if k == maxiter:
warnings.warn("Maximum number of iterations reached.")
return x[:k+1]About the code
The output of Function 4.5.2 float before the iteration starts.
Example 4.5.3 (Convergence of Newton’s method for systems)
Example 4.5.3 A system of nonlinear equations is defined by its residual and Jacobian.
Tip
Be careful when coding a Jacobian all in one statement. Spaces separate columns, so x[3]-1 is not the same as x[3] - 1.
function func(x)
[exp(x[2] - x[1]) - 2,
x[1] * x[2] + x[3],
x[2] * x[3] + x[1]^2 - x[2]
]
end;
function jac(x)
[
-exp(x[2] - x[1]) exp(x[2] - x[1]) 0
x[2] x[1] 1
2*x[1] x[3]-1 x[2]
]
end;
We will use a BigFloat starting value, and commensurately small stopping tolerances, in order to get a sequence long enough to measure convergence.
x₁ = BigFloat.([0, 0, 0])
ϵ = eps(BigFloat)
x = FNC.newtonsys(func, jac, x₁, xtol=ϵ, ftol=ϵ);
Let’s compute the residual of the last result in order to check the quality.
r = x[end]
@show residual = norm(func(r));residual = norm(func(r)) = 0.0
We take the sequence of norms of errors, applying the log so that we can look at the exponents.
logerr = [Float64(log(norm(r - x[k]))) for k in 1:length(x)-1]
ratios = [NaN; [logerr[i+1] / logerr[i] for i in 1:length(logerr)-1]]
@pt :header=["iteration", "log error", "ratio"] [eachindex(logerr) logerr ratios]The ratio is neatly converging toward 2, which is expected for quadratic convergence.
Example 4.5.3 A system of nonlinear equations is defined by its residual and Jacobian.
Tip
This function needs to be defined within a script file or in a file of its own with the .m extension.
function [f, J] = nlsystem(x)
f = zeros(3, 1); % ensure a column vector output
f(1) = exp(x(2) - x(1)) - 2;
f(2) = x(1) * x(2) + x(3);
f(3) = x(2) * x(3) + x(1)^2 - x(2);
J(1, :) = [-exp(x(2) - x(1)), exp(x(2) - x(1)), 0];
J(2, :) = [x(2), x(1), 1];
J(3, :) = [2 * x(1), x(3)-1, x(2)];
endSince our system function is defined in an external file here, we need to use @ in order to reference it as a function argument.
nlsystem = @f45_nlsystem;
x1 = [0; 0; 0]; % column vector!
x = newtonsys(nlsystem, x1);
num_iter = size(x, 2)Unrecognized function or variable 'newtonsys'. Let’s compute the residual of the last result in order to check the quality.
r = x(:, end)
back_err = norm(nlsystem(r))We take the sequence errors in the first component of the solution, applying the log so that we can look at the exponents.
log10( abs(x(1, 1:end-1) - r(1)) )'
ans =
0x1 empty double column vector
This sequence looks to be nearly doubling at each iteration, which is a good sign of quadratic convergence.
Example 4.5.3 A system of nonlinear equations is defined by its residual and Jacobian.
def func(x):
return array([
exp(x[1] - x[0]) - 2,
x[0] * x[1] + x[2],
x[1] * x[2] + x[0]**2 - x[1]
])
def jac(x):
return array([
[-exp(x[1] - x[0]), exp(x[1] - x[0]), 0],
[x[1], x[0], 1],
[2 * x[0], x[2] - 1, x[1]],
])
Our initial guess at a root is the origin.
x1 = zeros(3)
x = FNC.newtonsys(func, jac, x1)
print(x)[[ 0.00000000e+00 0.00000000e+00 0.00000000e+00]
[-1.00000000e+00 -1.20292797e-16 0.00000000e+00]
[-5.78586294e-01 1.57172588e-01 1.57172588e-01]
[-4.63138615e-01 2.30903685e-01 1.15452497e-01]
[-4.58026868e-01 2.35120714e-01 1.07713160e-01]
[-4.58033281e-01 2.35113900e-01 1.07689991e-01]
[-4.58033281e-01 2.35113900e-01 1.07689991e-01]]
The output has one row per iteration, so the last row contains the final Newton estimate. Let’s compute its residual.
r = x[-1]
f = func(r)
print("final residual:", f)final residual: [0.00000000e+00 1.38777878e-17 0.00000000e+00]
Let’s check the convergence rate:
logerr = [log(norm(x[k] - r)) for k in range(x.shape[0] - 1)]
for k in range(len(logerr) - 1):
print(logerr[k+1] / logerr[k])0.7937993447128695
3.6959808854483063
2.4326597889977193
2.3110932374366495
2.132541114533629
The ratio is apparently converging toward 2, as expected for quadratic convergence.
4.5.4 Exercises ¶ Two curves in the ( u , v ) (u,v) ( u , v ) u log u + v log v = − 0.3 u\log u + v \log v = -0.3 u log u + v log v = − 0.3 u 4 + v 2 = 1 u^4 + v^2 = 1 u 4 + v 2 = 1
(a) ✍ Write the intersection of these curves in the form f ( x ) = 0 \mathbf{f}(\mathbf{x}) = \boldsymbol{0} f ( x ) = 0 f \mathbf{f} f x \mathbf{x} x
(b) ✍ Find the Jacobian matrix of f \mathbf{f} f
(c) ⌨ Use Function 4.5.2 u = 1 u=1 u = 1 v = 0.1 v=0.1 v = 0.1
(d) ⌨ Use Function 4.5.2 u = 0.1 u=0.1 u = 0.1 v = 1 v=1 v = 1
Two elliptical orbits ( x 1 ( t ) , y 1 ( t ) ) (x_1(t),y_1(t)) ( x 1 ( t ) , y 1 ( t )) ( x 2 ( t ) , y 2 ( t ) ) (x_2(t),y_2(t)) ( x 2 ( t ) , y 2 ( t ))
[ x 1 ( t ) y 1 ( t ) ] = [ − 5 + 10 cos ( t ) 6 sin ( t ) ] , [ x 2 ( t ) y 2 ( t ) ] = [ 8 cos ( t ) 3 + 12 sin ( t ) ] , \begin{bmatrix}
x_1(t) \\ y_1(t)
\end{bmatrix}
=
\begin{bmatrix}
-5+10\cos(t) \\ 6\sin(t)
\end{bmatrix}, \qquad
\begin{bmatrix}
x_2(t)\\y_2(t)
\end{bmatrix} =
\begin{bmatrix}
8\cos(t) \\ 3+12\sin(t)
\end{bmatrix}, [ x 1 ( t ) y 1 ( t ) ] = [ − 5 + 10 cos ( t ) 6 sin ( t ) ] , [ x 2 ( t ) y 2 ( t ) ] = [ 8 cos ( t ) 3 + 12 sin ( t ) ] , where t t t
(a) ⌨ Make a plot of the two orbits with the following code:
x1(t) = -5+10*cos(t); y1(t) = 6*sin(t);
plot(x1,y1,0,2pi,aspect_ratio=1,legend=false)
x2(t) = 8*cos(t); y2(t) = 3 + 12*sin(t);
plot!(x2,y2,0,2pi)(b) ✍ Write out a 2 × 2 2\times 2 2 × 2
(c) ✍ Write out the Jacobian matrix of this nonlinear system.
(d) ⌨ Use Function 4.5.2
⌨ Suppose one wants to find the points on the ellipsoid x 2 / 25 + y 2 / 16 + z 2 / 9 = 1 x^2/25 + y^2/16 + z^2/9 = 1 x 2 /25 + y 2 /16 + z 2 /9 = 1 ( 5 , 4 , 3 ) (5,4,3) ( 5 , 4 , 3 )
x − 5 = λ x 25 , y − 4 = λ y 16 , z − 3 = λ z 9 , 1 = 1 25 x 2 + 1 16 y 2 + 1 9 z 2 \begin{split}
x-5 &= \frac{\lambda x}{25}, \\[1mm]
y-4 &= \frac{\lambda y}{16}, \\[1mm]
z-3 &= \frac{\lambda z}{9}, \\[1mm]
1 &= \frac{1}{25}x^2 + \frac{1}{16}y^2 + \frac{1}{9}z^2
\end{split} x − 5 y − 4 z − 3 1 = 25 λ x , = 16 λ y , = 9 λ z , = 25 1 x 2 + 16 1 y 2 + 9 1 z 2 for an unknown value of λ.
(a) Write out this system in the form f ( u ) = 0 \mathbf{f}(\mathbf{u}) = \boldsymbol{0} f ( u ) = 0
(b) Write out the Jacobian matrix of this system.
(c) Use Function 4.5.2 ( 5 , 4 , 3 ) (5,4,3) ( 5 , 4 , 3 )
⌨ Any three noncollinear points in the plane determine a unique circle. Suppose the points are given as ( x i , y i ) (x_i,y_i) ( x i , y i ) i = 1 , 2 , 3 i=1,2,3 i = 1 , 2 , 3 ( a , b ) (a,b) ( a , b ) r r r
f i ( a , b , r ) = ( a − x i ) 2 + ( b − y i ) 2 − r 2 f_i(a,b,r) = (a-x_i)^2 + (b-y_i)^2 - r^2 f i ( a , b , r ) = ( a − x i ) 2 + ( b − y i ) 2 − r 2 should be made zero for all i = 1 , 2 , 3 i=1,2,3 i = 1 , 2 , 3 f ( v ) = 0 \mathbf{f}(\mathbf{v})=\boldsymbol{0} f ( v ) = 0 v = [ a , b , r ] \mathbf{v}=[a,b,r] v = [ a , b , r ]
Use Function 4.5.2 ( − 5 , 0 ) (-5,0) ( − 5 , 0 ) ( 1 , − 3 ) (1,-3) ( 1 , − 3 ) ( 4 , 2 ) (4,2) ( 4 , 2 )