Usage of newtonsys

As before, a system of nonlinear equations is defined by its residual and Jacobian.

using FundamentalsNumericalComputation
function nlfun(x)
    f = zeros(3)  
    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];
    return f
end
   
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)

Our initial root estimate is the origin. The output is a vector of vector root estimates.

x1 = [0,0,0]

x = FNC.newtonsys(nlfun,nljac,x1)
7-element Array{Array{Float64,1},1}:
 [0.0, 0.0, 0.0]
 [-1.0, 0.0, 0.0]
 [-0.578586294114295, 0.1571725882285898, 0.15717258822858976]
 [-0.4631386148886986, 0.23090368503021674, 0.11545249687009654]
 [-0.4580268675321411, 0.23512071352763037, 0.10771316029325663]
 [-0.45803328074878114, 0.2351138998112444, 0.10768999092388234]
 [-0.45803328064126886, 0.23511389991867648, 0.10768999090411435]

We’ll compute the residual of the last result in order to check the quality.

r = x[end]
@show residual = nlfun(r);
residual = nlfun(r) = [0.0, 1.3877787807814457e-17, 0.0]

We take the sequence of norms of errors, applying the log so that we can look at the exponents.

[ log(norm(x[k]-r)) for k in 1:length(x)-1 ]
6-element Array{Float64,1}:
  -0.6424646089199801
  -0.5099879855618902
  -1.8849058464450312
  -4.585334658693529
 -10.597135921091917
 -22.59882803634857

The exponents approximately double, as is expected of quadratic convergence.