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.

Fundamentals of Numerical Computation

The Newton idea for systems¶