Building numerical jacobian matrix with complex functions (Aircraft Trim Code)

I am trying to build a 10X10 Jacobian matrix numerically with complex functions that call many other functions. Instead of doing each of the 100 partial derivatives seperately I'd like to design a code that will calculate each and build the matrix.

The code below works when the functions just use one variable, but it doesn't carry the other variables through the function and I get a "not enough input arguments error".

x = [Omegan;Omegas;Omegae;Omegaw;theta;phi;CTn;CTs;CTe;CTw];

F = [Fx(Omegan,Omegas,Omegae,Omegaw,theta,phi,CTn,CTs,CTe,CTw,u,v,w);

Fy(Omegan,Omegas,Omegae,Omegaw,theta,phi,CTn,CTs,CTe,CTw,u,v,w);

Fz(Omegan,Omegas,Omegae,Omegaw,theta,phi,CTn,CTs,CTe,CTw,u,v,w);

Mxt(Omegan,Omegas,Omegae,Omegaw,theta,phi,CTn,CTs,CTe,CTw,u,v,w);

Myt(Omegan,Omegas,Omegae,Omegaw,theta,phi,CTn,CTs,CTe,CTw,u,v,w);

Mzt(Omegan,Omegas,Omegae,Omegaw,theta,phi,CTn,CTs,CTe,CTw,u,v,w);

CT(Omegan,CTn,rho,A,R,u,v,w,Xin,Yin,kn,theta);

CT(Omegas,CTs,rho,A,R,u,v,w,Xis,Yis,ks,theta);

CT(Omegae,CTe,rho,A,R,u,v,w,Xie,Yie,ke,theta);

CT(Omegaw,CTw,rho,A,R,u,v,w,Xiw,Yiw,kw,theta)];

for k = 1:numel(x)
    for j = 1:numel(x)
        Jacobian(k,j) = Derivative(eval(['F' num2str(k)]),x,x(j),j);
    end
end

where,

function J = Derivative(F,x,xi,i)

e = 0.0001;
x1 = x;
x2 = x;
x1(i) = xi;
x2(i) = xi+e;
J = (F(x2)-F(x1))./e;

end

Is there a way for me to build the numerical jacobian with a script rather than calculating the numerical derivative for each with the derivative formula?

For background, I am using newtons method to solve for [Omegan;Omegas;Omegae;Omegaw;theta;phi;CTn;CTs;CTe;CTw] at various speeds u , v , w

Answer 1

[Gobiha Duraisamy]

The jacobian function could be used to compute the Jacobian matrix. Use syms to create symbolic state variables,

syms Omegan Omegas Omegae Omegaw theta phi CTn CTs CTe CTw

Define vector functions (Fx,Fy,Fz,Mxt,Myt,Mzt,CT1,CT2,CT3,CT4) using the symbolic variables. Now compute the Jacobian matrix using the following syntax:

j = jacobian([Fx,Fy,Fz,Mxt,Myt,Mzt,CT1,CT2,CT3,CT4],[Omegan,Omegas,Omegae,Omegaw,theta,phi,CTn,CTs,CTe,CTw])

Use subs to substitute symbolic variables with a value,

s = subs(j, [Omegan,Omegas,Omegae,Omegaw,theta,phi,CTn,CTs,CTe,CTw], [1,1,2,2,2,3,4,5,5,5]);
sNum = double(s)