Using levenberg¶
To solve a nonlinear system, we need to code only the function defining the system (residual vector), and not its Jacobian.
using FundamentalsNumericalComputation
function nlsystem(x)
return [
exp(x[2]-x[1]) - 2,
x[1]*x[2] + x[3],
x[2]*x[3] + x[1]^2 - x[2]
]
end
nlsystem (generic function with 1 method)
In all other respects usage is the same as for the newtonsys function.
x1 = [0,0,0]
x = FNC.levenberg(nlsystem,x1)
12-element Array{Array{Float64,1},1}:
[0.0, 0.0, 0.0]
[-0.08396946555420857, 0.07633587776621241, 0.0]
[-0.4220507601545548, 0.21991260524972361, 0.012997569767382334]
[-0.48610710902101084, 0.21389682964709655, 0.0977187252071678]
[-0.45628390779216876, 0.2421104774129355, 0.10100440266654645]
[-0.45563883362225116, 0.23470443539126568, 0.108546657255663]
[-0.4583961451217119, 0.23530956863847285, 0.10739828072313906]
[-0.4580434038157383, 0.23512124061139036, 0.10768079583181711]
[-0.4580333258441392, 0.23511390840122376, 0.1076899804953973]
[-0.45803327880719313, 0.23511389867393445, 0.10768999250671268]
[-0.45803328056019965, 0.23511389986307896, 0.107689990975689]
[-0.458033280641234, 0.2351138999186528, 0.10768999090414473]
It’s always a good idea to check the accuracy of the root, by measuring the residual (backward error).
r = x[end]
@show backward_err = norm(nlsystem(r));
backward_err = norm(nlsystem(r)) = 1.2708308198538738e-13
Looking at the convergence of the first component, we find a subquadratic convergence rate, just as with the secant method.
[ log(norm(x[k]-r)) for k in 1:length(x)-1 ]
11-element Array{Float64,1}:
-0.642464608920046
-0.8665622983300257
-2.278548436446176
-3.3083905468408217
-4.621899500113314
-5.961582353290289
-7.590899766347262
-11.073290988317947
-16.86982474417951
-19.71706564311687
-22.830722547915087