I'm trying to write a function that accepts a function vector of any size and return the numerical Jacobian matrix.

The answer should be x=1.2812, y=2.5349 and z=-0.5022 but once I put the loop, I'm getting these values x=0.9794, y=2.9794 and z=-0.0206

This is my function script:

function J = jacobian_numeric(f,x,h)

[m,n] = size(x);

for ii = 1:n
    hv = zeros(1,n);
    hv(1,ii) = h;
    J(:,ii) = (f(x+hv) - f(x-hv))./(2*h);
end
end

And this is the m-file that calls the function:

clear; clc;
f = @(x) [(x(1).^2 - x(2) - x(1).*x(3) + 0.25) ; (x(1).*x(2) + x(3).^2 - 3.5) ; (sin(pi.*x(1)) - x(2).*x(3) - 0.5)];

x0 = [1 ; 3 ; 0];
h = 0.00001;
tol = 1e-5;
err = 1 + tol;
iter = 0;
while err > tol
    J = jacobian_numeric(f,x0,h);
    xold = x0;
    x0 = x0 - J \ f(x0);
    iter = iter + 1;
    err = max(abs((x0 - xold)./(x0)));
end

x = x0(1,1)
y = x0(2,1)
z = x0(3,1)

I tried many things but nothing worked, so I am hoping on getting help from someone here. Thank you in advance.

[Torsten]

In jacobian_numeric, how should MATLAB know that the first dimension of J has to be n ? Initialize J as J = zeros(n,n) !

Best wishes

Torsten.

[bhador]

I don't get what you mean, and I tried adding J = zeros(n,n) but did not work.

[Torsten]

function J = jacobian_numeric(f,x,h)

n = numel(x);
J = zeros(n,n);

for ii = 1:n
    hv = zeros(1,n);
    hv(1,ii) = h;
    J(:,ii) = (f(x+hv) - f(x-hv))./(2*h);
end
end

[bhador]

The script is not running, it says error. And the answer I'm getting is NaN

[Torsten]

function J = jacobian_numeric(f,x,h)

n = numel(x);
J = zeros(n,n);

for ii = 1:n
    hv = zeros(n,1);
    hv(ii,1) = h;
    J(:,ii) = (f(x+hv) - f(x-hv))./(2*h);
end
end

Wouldn't have happened without implicit expansion ...

Best wishes

Torsten.