How do I evaluate a symbolic Jacobian at multiple values of the same variable in the same equation?

Hi @Andrew Relyea,

You mentioned, “I need the function evaluated at L, l, a, r, p, y and Li, li, ai, ra, pa, ya but since the later aren't different states I can't just do

statevect = [r, p, y, L, l, a, ra, pa, ya, Li, li, ai]

because it would differentiate the sets of state values as separate states, and

statevect = [r, p, y, L, l, a]

won't work because any term with the second set of terms won't be differentiated. Any help would be greatly appreciated.”

Please see my response to your comments below.

You are working with a complex system where you have two different time steps represented by two sets of state variables and your goal is to compute the Jacobian of a derived function (feat_pos_g) that combines these states in a specific transformation. However, you are facing a challenge since the symbolic differentiation needs to consider both sets of variables simultaneously. After doing analysis of your code and reviewing the jacobian documentation provided at the link below,

https://www.mathworks.com/help/symbolic/sym.jacobian.html

here is a revised version that correctly constructs the state vector and calculates the Jacobian:

syms L l a r p y Rn Re Li li ai ra pa ya Rei az el azi eli real

% Constants
E2 = 6.694379990141317e-3; % Eccentricity constant
kappa = 1 - E2;

% Transformation from Lla to NED Cartesian
Ce2n = [-sin(L)*cos(l), -sin(L)*sin(l), cos(L); 
      -sin(l), cos(l), 0; 
      -cos(L)*cos(l), -cos(L)*sin(l), -sin(L)];

p0 = Ce2n * [(Re + a) * cos(L) * cos(l); 
            (Re + a) * cos(L) * sin(l); 
            (kappa * Re + a) * sin(L)];

pj = Ce2n * [(Rei + ai) * cos(Li) * cos(li); 
            (Rei + ai) * cos(Li) * sin(li); 
            (kappa * Rei + ai) * sin(Li)];

% Create the unit LOS vector
cr = cos(r); sr = sin(r);
cp = cos(p); sp = sin(p);
cy = cos(y); sy = sin(y);

Cb2n0 = [cp*cy, cp*sy, -sp; 
       -cr*sy + sr*sp*cy, cr*cy + sr*sp*sy, sr*cp; 
       sr*sy + cr*sp*cy, -sr*cy + cr*sp*sy, cr*cp]';

cr2 = cos(ra); sr2 = sin(ra);
cp2 = cos(pa); sp2 = sin(pa);
cy2 = cos(ya); sy2 = sin(ya);

Cb2nj = [cp2*cy2, cp2*sy2, -sp2; 
        -cr2*sy2 + sr2*sp2*cy2, cr2*cy2 + sr2*sp2*sy2, sr2*cp2; 
        sr2*sy2 + cr2*sp2*cy2, -sr2*cy2 + cr2*sp2*sy2, cr2*cp2]';

u0 = Cb2n0 * [cos(az)*cos(el);
             sin(az)*cos(el);
             sin(el)];
uj = Cb2nj * [cos(azi)*cos(eli);
             sin(azi)*cos(eli);
             sin(eli)];

% Trilaterate to initialize feature position
p1 = dot(cross(pj - p0, uj), cross(u0, uj)) / (norm(cross(u0, uj))^2);
p2 = dot(cross(p0 - pj, u0), cross(uj, u0)) / (norm(cross(uj, u0))^2);

feat_pos_n = 1/2 * (p0 + pj + p1 .* u0 + p2 .* uj);

feat_pos_g = [L; l; a] + [1/(Rn + a), 0, 0; 
                        0, 1/((Re + a) * cos(L)), 0; 
                        0, 0, -1] * (feat_pos_n - p0);

V% Create a combined state vector for differentiation
statevect = [L l a r p y Li li ai ra pa ya];

% Calculate the Jacobian of feat_pos_g with respect to the state vector
J_feat_pos_g = jacobian(feat_pos_g, statevect);

% Display the Jacobian
disp(J_feat_pos_g);

Explanation of modified code

Symbolic Variables: define all necessary symbolic variables at once.

Coordinate Transformations: The code maintains your transformations from ECEF to NED Cartesian coordinates.

Unit Vectors: The unit line-of-sight vectors are created using rotation matrices as per your original structure.

Feature Position Calculation: The trilateration method remains unchanged but is now part of a clean computation pipeline.

State Vector Construction: A single combined statevect is constructed from both sets of states ([L l a r p y Li li ai ra pa ya]).This allows MATLAB's jacobian function to differentiate correctly across all required variables.

Jacobian Calculation: Finally, computing the Jacobian using MATLAB's built-in jacobian function and display it.

Please see attached.

Make sure that you validate the outputs from each transformation step for consistency. Now, if performance becomes an issue with more complex systems or larger symbolic expressions, consider numerical approximations or simplifying assumptions where applicable. This approach should meet your requirements and help you effectively navigate through your symbolic differentiation task. If you have any further questions, please let me know.

Hi @Andrew Relyea,

You mentioned, “I have attached a screen capture to show the form I need the output matrix to take. As you can see, I am using the output to augement another matrix. In the original there is not, for example, a state Li but there is a L. Creating additional columns misaligns the output.”

Please see my response to your comments below.

After reviewing your attached a screen capture to show the form you need the output matrix to take and reading your comments, please follow my suggestions mentioned below.

Instead of separating your states into two different groups, create a unified state vector that includes both sets of variables. For example:

statevect = [r; p; y; L; l; a; ra; pa; ya; Li; li; ai];

Make sure that your function feat_pos_g is defined in such a way that it takes this unified vector as input. This will allow you to evaluate all states simultaneously:

   function feat_pos_g = calculate_feat_pos_g(statevect)
       % Extract states
       r = statevect(1);
       p = statevect(2);
       y = statevect(3);
       L = statevect(4);
       l = statevect(5);
       a = statevect(6);
       ra = statevect(7);
       pa = statevect(8);
       ya = statevect(9);
       Li = statevect(10);
       li = statevect(11);
       ai = statevect(12);

       % [Insert calculations for feat_pos_n and feat_pos_g here]
   end

Then, use automatic differentiation tools available in MATLAB or numerical differentiation methods to compute the Jacobian. Here’s an example using finite differences:

   function J = calculate_jacobian(statevect)
       epsilon = 1e-6; % Small perturbation for finite difference
       n = length(statevect);
       m = 3; % Assuming feat_pos_g has 3 outputs

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

       for i = 1:n
           perturb = zeros(n, 1);
           perturb(i) = epsilon;

           f_plus = calculate_feat_pos_g(statevect + perturb);
           f_minus = calculate_feat_pos_g(statevect - perturb);

           J(:, i) = (f_plus - f_minus) / (2 * epsilon);
       end
   end

When using this Jacobian in subsequent calculations (like augmenting another matrix), you have to make sure that it matches the expected dimensions as indicated in your example matrices. For instance, if you expect a specific structure like Paug, make sure JGp and JGz are appropriately sized.

After implementing these changes, validate your results by comparing against known cases or using test scenarios where the expected outcomes are established.

Hope this helps.