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.