Untitled

We repeat a previous example. We need to code only the function defining the system, and not its Jacobian. The function is found at the end of this file.

In all other respects usage is the same as for the newtonsys function.

x1 = [0;0;0];   
x = levenberg(@nlsystem,x1);

It's always a good idea to check the accuracy of the root, by measuring the residual (backward error).

r = x(:,end)
r = 
   -0.4580
    0.2351
    0.1077
backward_err = norm(nlsystem(r))
backward_err = 1.2708e-13

Looking at the convergence of the first component, we find a subquadratic convergence rate, just as with the secant method.

log10( abs(x(1,1:end-1)-r(1)) )'
ans = 
   -0.3391
   -0.4271
   -1.4439
   -1.5517
   -2.7571
   -2.6208
   -3.4403
   -4.9947
   -7.3448
   -8.7366

function f = nlsystem(x)
    f = zeros(3,1);   % ensure a column vector output
    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);
end