# 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 = exp(x-x) - 2;
f = x*x + x;
f = x*x + x^2 - x;
return f
end

function nljac(x)
J = zeros(3,3)
J[1,:] = [-exp(x-x),exp(x-x), 0]
J[2,:] = [x, x, 1]
J[3,:] = [2*x, x-1, x]
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.