Unfortunately, your update doesn't convince me you did the correct job, based on your script, your x variable is already the answer and residing in the workspace, the second time when your own scripts are being executed, only one iteration has been made, not the whole Newton's method.
If you are allowed to use Symbolic Toolbox to compute the Jacobian, there is a more universal way of doing instead of just carrying out a problem-oriented script. Here is a modified version to match your notation of an old implementation of mine for Newton's method, and this could be easily vectorized for a multi-dimensional nonlinear equation system using varargin input, and do a string size check on the inline function you passed to the following function.
function [s J] = newton(f,p0,tol,MaxIter)
format long
x = sym('x'); y = sym('y');
F = f([x,y]);
% Compute the Jacobian matrix symbolically
J = jacobian(F);
invJ = inv(J);
s = zeros(MaxIter,2);
s(1,:) = p0;
dsnorm = inf;
iter = 1;
while dsnorm>tol && iter<MaxIter
ds = -subs(invJ,[x y],s(iter,:))*f(s(iter,:));
s(iter+1,:) = s(iter,:) + ds';
dsnorm = norm(ds,inf);
iter = iter+1;
end
s = s(1:iter,:);
end
This function accepts input f as an inline function, p0 is the initial guess, tol is the tolerance and MaxIter is the maximum iteration allowed, if the vectorized version of your input is:
>>f=inline('[p(:,1).^2+p(:,2).^2-2.12; ...
p(:,2).^2-p(:,1).^2.*p(:,2)-0.04]','p');
Then
>>s = newton(f,[1 1],0.4e-6,10)
would return:
s =
1.000000000000000 1.000000000000000
1.006666666666667 1.053333333333333
1.006852189012140 1.051784722315402
1.006852720492674 1.051783057117558
1.006852720493521 1.051783057115295
As you can see the last output is different from yours, because basically what you had was fsolve output, however fsolve use trust-region-dogleg algorithm by default, not Newton's method, if your output coincides with fsolve's result, then your implementation should have a problem.
As
$$
f(x) = f(x_0)+\frac{\partial f(x_0)}{\partial x}(x-x_0) + \frac 12(x-x_0)^{\dagger}\frac{\partial^2f(x_0)}{\partial^2x}(x-x_0) + O(|x-x_0|^3)
$$
deriving again we have
$$
\nabla f(x)\approx\nabla f(x_0) + \frac{\partial^2f(x_0)}{\partial^2x}(x-x_0)\approx 0
$$
because $f(x) = g_1(x)^2+g_2(x)^2+g_3(x)^2$ so we get
$$
x_k = x_{k-1}-H_{k-1}^{-1}\nabla f(x_{k-1})
$$
with
$$
\nabla f(x) = \frac{\partial f(x_0)}{\partial x} = J\\
H = \frac{\partial^2f(x)}{\partial^2x}\approx J^{\dagger}\cdot J
$$
NOTE
Calling $G(x) = \{g_1(x),g_2(x),\cdots,g_m(x)\}^{\dagger}$ and $f(x) = \frac 12G^{\dagger}(x)\cdot G(x)$ we have
$$
\nabla f(x) = J(x)^{\dagger}\cdot G(x)
$$
and
$$
\nabla^2 f(x) = H(x) = \nabla J(x)^{\dagger}\cdot G(x) + J(x)^{\dagger}\cdot \nabla G(x) = J(x)^{\dagger}\cdot J(x) +\sum_kg_k(x)\nabla^2 g_k(x)
$$
So we can evaluate the approximation
$$
H = \frac{\partial^2f(x)}{\partial^2x}\approx J(x)^{\dagger}\cdot J(x)
$$
Attached a MATHEMATICA script implementing the Gauss-Newton iterative procedure
X = {x, y, z};
f[x_, y_, z_] := {{2 x - y z - 1}, {1 - x + y - Exp[x - z]}, {-x - 2 y + 3 z}}
g[x_, y_, z_] := D[f[x, y, z], {X}]
x0 = RandomReal[{-2, 2}];
y0 = RandomReal[{-2, 2}];
z0 = RandomReal[{-2, 2}];
X0 = {x0, y0, z0};
f0 = Flatten[f[x0, y0, z0], 1];
n = 20;
error = 0.0001;
path = {{0, x0, y0, z0, f0.f0}};
For[k = 1, k <= n, k++,
Jf = Flatten[Transpose[g[x, y, z] /. Thread[X -> X0]], 1];
JR = Transpose[Transpose[f[x, y, z]].Jf] /. Thread[X -> X0];
Hf = Transpose[Jf].Jf;
p = LinearSolve[Hf, -JR];
If[Max[Abs[p]] < error, Break[]];
X0 = Flatten[X0 + p];
f0 = Flatten[f[x, y, x], 1] /. Thread[X -> X0];
AppendTo[path, {k, x, y, z, f0.f0} /. Thread[X -> X0]]
]
with a typical outcome
$$
\left[
\begin{array}{ccccc}
k & x & y & z & f^{\dagger}\cdot f\\
0 & -0.18083 & -0.360527 & -1.54181 & 27.0262 \\
1 & 0.419264 & -0.249572 & -0.0266268 & 2.23994 \\
2 & 0.441999 & 0.376084 & 0.398056 & 0.101378 \\
3 & 0.693323 & 0.692737 & 0.692933 & 0.00877097 \\
4 & 0.846417 & 0.846417 & 0.846417 & 0.000556373 \\
5 & 0.923209 & 0.923209 & 0.923209 & 0.0000347735 \\
6 & 0.961604 & 0.961604 & 0.961604 & \text{2.1733}\cdot 10^{-6}
\\
7 & 0.980802 & 0.980802 & 0.980802 & \text{1.3583}\cdot 10^{-7} \\
8 & 0.990401 & 0.990401 & 0.990401 & \text{8.4896}\cdot 10^{-9} \\
9 & 0.995201 & 0.995201 & 0.995201 & \text{5.3060}\cdot 10^{-10}
\\
10 & 0.9976 & 0.9976 & 0.9976 & \text{3.3162}\cdot 10^{-11} \\
11 & 0.9988 & 0.9988 & 0.9988 & \text{2.0726}\cdot 10^{-12} \\
12 & 0.9994 & 0.9994 & 0.9994 & \text{1.2954}\cdot 10^{-13} \\
13 & 0.9997 & 0.9997 & 0.9997 & \text{8.0963}\cdot 10^{-15} \\
14 & 0.99985 & 0.99985 & 0.99985 & \text{5.0603}\cdot 10^{-16} \\
\end{array}
\right]
$$
Attached a script in octave
function [X0,k] = GaussNewton(x0,y0,z0,n,error)
format long
%n = 20;
%error = 0.0001;
%x0 = 0;
%y0 = 0;
%z0 = 0;
X0 = [x0;y0;z0];
for k = 1:n
x0 = X0(1);
y0 = X0(2);
z0 = X0(3);
Jf = g(x0,y0,z0);
JR = (f(x0,y0,z0)'*Jf)';
Hf = Jf'*Jf;
p =-inv(Hf)*JR;
if max(abs(p)) < error
break;
end
X0 = X0 + p;
end
end
function fx = f(x,y,z)
fx = [2.*x-y*z-1.;1.-x+y-exp(x-z);-x-2.*y+3.*z];
end
function dfx = g(x,y,z)
dfx = [2.,-z,-y;-exp(x-z)-1.,1.,exp(x-z);-1.,-2.,3.];
end
Best Answer
So you want to write a subroutine that returns the jacobian matrix at point $x\in \mathbb{R}^2$ of a function handle
@fin matlab. Also you wanna do this numerically, not symbolically.Here is a simple finite difference subroutine that does the job for
fbeing a function handle:The input
fis a function handle,xis any point. The outputJis the jacobian matrix atx. Example: $f = (x^2+y^2, x^2-y^2)$:For the
Symbolic Toolboxapproach, please refer to my answer here in this question: Solving a system with Newton's method in matlab? . A more object-oriented approach would be like you mentioned, allowing user to input any type of function: inline function, symbolic, function handle, etc. Converting symbolic function to function handle can be done bymatlabFunctionin MATLAB. And inline functions arestring, hence easy to be converted to a symbolic function.