The Newton idea for systems

Let us use Newton’s method on the system defined by the function

using LinearAlgebra
function nlvalue(x)
    return [ exp(x[2]-x[1]) - 2,
          x[1]*x[2] + x[3],
          x[2]*x[3] + x[1]^2 - x[2]
        ]
end
nlvalue (generic function with 1 method)

Here is a function that computes the Jacobian matrix.

function nljac(x)
    J = zeros(3,3)
    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]]
    return J
end
nljac (generic function with 1 method)

(These functions could be written as separate files, or embedded within another function as we have done here.) Our initial guess at a root is the origin.

x = [0.,0.,0.]
3-element Array{Float64,1}:
 0.0
 0.0
 0.0

We need the value (residual) of the nonlinear function, and its Jacobian, at this value for \(\mathbf{x}\).

f = nlvalue(x)
3-element Array{Float64,1}:
 -1.0
  0.0
  0.0
J = nljac(x)
3×3 Array{Float64,2}:
 -1.0   1.0  0.0
  0.0   0.0  1.0
  0.0  -1.0  0.0

We solve for the Newton step and compute the new estimate.

s = -(J\f)
newx = x + s
3-element Array{Float64,1}:
 -1.0
  0.0
  0.0

Here is the new residual.

nlvalue(newx)
3-element Array{Float64,1}:
 0.7182818284590451
 0.0
 1.0

We don’t seem to be especially close to a root. Let’s iterate a few more times. To do that, we start keeping the root estimates as elements of a vector (i.e., a vector of vectors). Same for the residuals.

x = [x,newx]
2-element Array{Array{Float64,1},1}:
 [0.0, 0.0, 0.0]
 [-1.0, 0.0, 0.0]
y = [f,nlvalue(newx)]
2-element Array{Array{Float64,1},1}:
 [-1.0, 0.0, 0.0]
 [0.7182818284590451, 0.0, 1.0]
for n = 2:7
    f = nlvalue(x[n])
    push!(y,f)
    s = -nljac(x[n])\f
    push!(x,x[n]+s)
end

We find the infinity norms of the residuals.

residual_norm = [ norm(f,Inf) for f in y ]   # max in dimension 1
8-element Array{Float64,1}:
 1.0
 1.0
 1.0
 0.20229273399879752
 0.010252098544050403
 2.1556384274173945e-5
 1.9899884518004285e-10
 1.3877787807814457e-17

We don’t know an exact answer, so we will take the last computed value as its surrogate.

r = x[end]
3-element Array{Float64,1}:
 -0.45803328064126886
  0.23511389991867648
  0.10768999090411434

The following will subtract r from every column of x.

error = [ norm(xk-r) for xk in x[1:end-1] ]
7-element Array{Float64,1}:
 0.525994454581891
 0.6005027934725727
 0.15184335534811869
 0.010200335585503865
 2.4987473439459857e-5
 1.5326877679492286e-10
 1.3877787807814457e-17

Now we take infinity norms and check for quadratic convergence.

ratios = [ log(error[k+1]) / log(error[k]) for k in 1:length(error)-1 ]
6-element Array{Float64,1}:
 0.7937993447128695
 3.695980885448303
 2.4326597889977144
 2.3110932374368462
 2.1325411123294264
 1.7176219080468489

For a brief time, we see ratios around 2, but then the limitation of double precision makes it impossible for the doubling to continue.