How can I loop over i and j, to compute the ith row and jth column of the Jacobian matrix, using the central difference method?

Hi there!

I'm practicing with this toy problem and want to loop over i and j to compute the ith row and jth column of the Jacobian partial derivatives matrix, using the central difference method. The for loop that's commented out works already (it closely follows an answer posted yesterday by Matt J.), so I'm just trying to find another way to do this problem. Ideally, I want to do something like:

J(i,j) = ( zdot(i) ( X(j) + h - zdot(i) ( X(j) - h ) ) / (2*h), but I get errors from Matlab about chain indexing not being allowed.

Most natural for me would be to differentiate across, not vertically, when computing the Jacobian, so I was wondering whether I could put this into code, too.

Thanks in advance,

%% Jacobian Practice
zdot         = @(z) myrhs(z);
InitialGuess = [1,1];            % pass a good initial guess to fsolve
X = fsolve(zdot,InitialGuess)    % a fixed point
h = .000001;                     % finite-differencing step-size
% delta = speye(2);
% J = zeros(2,2);
% for i = 1:2
%     J(:,i) = ( zdot( X + h*delta(:,i) ) - zdot( X - h*delta(:,i) ) ) / (2*h);
% end
for i = 1:2

    for j = 1:2 

        J(i,j) = ( zdot( X(j) + h ) - zdot( X(j) - h) ) / (2*h);

    end

end
%% write a system of two equations in two unknowns
function zdot = myrhs(z)
x = z(1);
y = z(2);
xdot = x + y^2 - pi; 
ydot = x + y;
zdot = [xdot; ydot];
end

[Torsten]

X(j) + h is a single value while your function expects a vector with two elements.

[Noob]

Hi Torsten!

Yes, I've fixed that error. I've tried a number of different things but can't seem to loop over i and j, which would be nice to do, in order to have a second way of computing the Jacobian. So far, Matt's answer seems best to go with. I'll keep trying :)

Thanks!

[Torsten]

Evaluating "myrhs" at X+h*delta(:,i) and X-h*delta(:,i) suffices to get the complete i-th column of the Jacobian. So there is no need to loop over j. Of course, you could do that (see below), but why ?

%% Jacobian Practice
zdot         = @(z) myrhs(z);
InitialGuess = [1,1];            % pass a good initial guess to fsolve
X = fsolve(zdot,InitialGuess)    % a fixed point

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

X = 1×2

-2.3416 2.3416

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

X = X(:);
h = .000001;                     % finite-differencing step-size
J = zeros(2,2);
e = zeros(2,1);
for i = 1:2
    e(i) = 1;
    fp = zdot( X + h*e );
    fm = zdot( X - h*e );
    for j = 1:2 
        J(j,i) = ( fp(j) - fm(j) ) / (2*h);
    end
    e(i) = 0;
end
J

J = 2×2

1.0000 4.6833 1.0000 1.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

J_ana = [1 2*X(2);1 1]

J_ana = 2×2

1.0000 4.6833 1.0000 1.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%% write a system of two equations in two unknowns
function zdot = myrhs(z)
x = z(1);
y = z(2);
xdot = x + y^2 - pi; 
ydot = x + y;
zdot = [xdot; ydot];
end

[Noob]

Hi Torsten!

Why do you set e(i) = 0 before ending the nested loop?

There's perhaps no good reason to loop over both i and j, other than for practice, and to see another way to do it.

I took a peek at your solution. Let me see if I can write it myself.

Thanks so much for your help!

[Noob]

Hi Torsten!

Here's my own code, inspired by your answer, and it seems to work, provided that I sent e(i) = 0 before the loop ends.

Do you know why I must set e(i) = 0 before ending the loop?

Now I'm happy: I know three ways to write a central difference method now: Matt's, yours, and mine.

Thanks so much!

for i = 1:2 
    e(i) = 1;
    FL = ( zdot( X + h * e ) - zdot( X - h * e ) ) / (2*h);
    J(:,i) = FL;
    e(i) = 0;
    
end

[Noob]

Hi Torsten!

You set e(i) = 0; otherwise e = (1,1) when the loop index is i = 2, which will give an error.

I understand it now.

Thanks so much!

[Noob]

Hi @Torsten

Feel free to move your comment to an answer, so I can accept it.

Goodnight and thanks again!

[Matt J]

Now I'm happy: I know three ways to write a central difference method now: Matt's, yours, and mine.

There's really no substantive difference between the method you've settled upon and the one I originally proposed. You've simply assigned delta(:,i) to a temporary variable e and the finite difference to a temporary variable FL.

J = nan(N);
delta = speye(N);

%Yours, with minor re-arrangement
for i = 1:2   
    e=delta(:,i);
    FL = ( zdot( X + h * e ) - zdot( X - h * e ) ) / (2*h);
    J(:,i) = FL;
end

%Original
for i = 1:2
     J(:,i) = ( zdot( X + h*delta(:,i) ) - zdot( X - h*delta(:,i) ) ) / (2*h);
end

Note that,with this re-arrangement, e will be a sparse vector, so evaluating h*e will be faster.

[Noob]

Hi Matt!

Yes, thanks for that!

I was wondering:

How exactly does the central difference method differentiate one variable while keeping the other variables constant? When I do it with formulas it’s clear. If there are no formulas for the functions, how does that happen? Sorry for the basic question.

Thanks!

[Torsten]

If there are no formulas for the functions, how does that happen?

How do you want to apply central differencing if there are no formulas for the functions ?

As you can see from the central differencing formula, you must be given values for zdot in at least 4 points around X along the 2 main axes to compute an approximation for the Jacobian (and in 12 points along the 6 main axes for your ODE application).