Computing a Jacobian Numerically using 5pt stencil approximation

Hi guys! I am trying to compute a jacobian numerically using a 5-pt stencil approximation. my array F contains two functions:

 x = [x1;x2];

 f = @(x1,x2) func1(x);
 g = @(x1,x2) func2(x);

 %Assigning functions to an array
 F(1,1) = {f};
 F(2,1) = {g};

Then, I am trying to compute my jacobian but am not getting proper results.

function [J,h] = jacob(F,x)

[n,m] = size(x);
h = zeros(n,1);

%Initialize Jacobian
J = zeros(n,n);

%Numerical computation of Jacobian using 5-pt stencil approximation
for i = 1:n
    for j = 1:m
        %If i == j, h takes the value of the step size
        if i == j
            h(i) = 1e-3;
        end         
        J(i,j) = (F{i,1}(x(j)+2*h(j)) + 8*F{i,1}(x(j)+h(j)) - 8*F{i,1}(x(j)-h(j)) + F{i,1}(x(j)-2*h(j)))/(12*h(j));
        h(j) = 0;
    end
end

end

[Matt J]

I am trying to compute my jacobian but am not getting proper results.

as demonstrated by...?

[Walter Roberson]

The question is how you know you are getting results that are not proper. What should we be looking at? If we were to make a change to your code in hopes of fixing the problem, then how would we know if we had succeeded ?

[Cassandra Athans]

Basically, the problem is that when using the function handle, the value of x = [x1,x2] is automatically fixed. In F = [f;g] where f = f(x1,x2) and g = f(x1,x2).

When I try evaluating F at x+h or x-h, the value of x will not change. I am therefore having trouble evaluating F at different values of x, or manipulating x itself.

Answer 1

[Walter Roberson]

 f = @(x1,x2) func1([x1,x2]);
 g = @(x1,x2) func2([x1,x2]);

However, in your code you invoke

F{i,1}(x(j)+2*h(j))

so you carefully defined a function handle to take two values, but you are passing in only one value.

[Cassandra Athans]

when invoking

F{i,1}(x(j)+2*h(j))

should I then pass in 2 values? It seems as if when I change what is in (...), the value of F{i,1} does not change. It remains constant @x = [x1;x2]

[Walter Roberson]

func1 = @(xy) xy(1).^2 - sin(xy(2));
func2 = @(xy) 2*xy(1) + cos(xy(2));

f = @(x1,x2) func1([x1,x2]);
g = @(x1,x2) func2([x1,x2]);

   F(1,1) = {f};
   F(2,1) = {g};

   F{1,1}(7,1/2)
   F{2,1}(7,1/2)

   syms x y
   F{1,1}(x, y)
   F{2,1}(x, y)

Remember, when you just look at F{1,1} you are just looking at the function handle, without having evaluated it. You need to pass arguments to get evaluation.