MATLAB - Newton's Method to solve a system of nonlinear equations

Question

I'm trying to write a function that uses Newton's Method to solve a system of nonlinear equations

[roots, count, resids, history] = Newtons(func,x0,tol)

with the following inputs and outputs.

I'm struggling with creating the function handle and using it in the main function as well as understanding and generating the resids and history outputs in particular. I would really appreciate some help with this function.

What I've done/tried so far: (I'm relatively new to Matlab and believe it is totally wrong and not worth looking at)

function [f, J] = jacobian(f, x)
n=length(x); 
fx=f(x); 
step=1e-10; 

for i=1:n
   xstep = x;
   xstep(i)=x(i)+step;
   J(:,i)=(f(xstep)-fx)/step;
end
end

function [roots, count, resids, history] = Newtons(func,x0,tol)
func = @jacobian;

if nargin == 2
    tol = 10e-10;
end

MAXIT = 50;
xx = x0;
N = 0;
while N < 50
    JJ = func(2);

end

Answer 1

In general one wants to avoid for-loops in Matlab and try and do as much of the code vectorized as possible as Matlab then can use some really fancy and optimised BLAS libraries for the computations.

The Jacobian is a Matrix which has one index for function input dimension number and the other for the partial derivative number.

If we decide to stack the returned function values along a third dimension, we can define functions like this :

 f = @(a)(cat(3,a(:,1)-a(:,2).^2,a(:,1).^3-a(:,2)));

This would be the function $(x,y)\to (x-y^2,x^3-y)$ with the input values in the different columns.

jacob_11 = reshape((f([1,1;1,1]+1e-3*eye(2))-f([1,1;1,1]))./(1e-3),[2,2])

This is the line to evaluate at the point $(1,1)$ . The step size here is the 1e-3

Just one line of code!

The output is: $$jacob_{11}=\left[\begin{array}{rr}1&3\\-2&-1\end{array}\right]$$ Which we can verify by hand since we used such a simple function.

---

Now the next step to solve the equation system $\cases{x-y^2=0\\x^3-y=0}$ we can do this by constructing the function $g = (x-y^2)^2+(x^3-y)^2$ and then minimizing it (why?)

Now since we have a scalar valued function our Jacobian becomes a gradient, which we can follow in a gradient descent. As we can see, we can manage to find the special solution $x=y=1$.

The next step towards a Newton method would be calculating a Hessian. That is matrix of partial second order derivatives to the objective function.