Jacobian of Massive Equations

4 views (last 30 days)
Ayden Clay
Ayden Clay on 1 Feb 2021
Edited: Ayden Clay on 1 Feb 2021
As the title states I am trying to find the Jacobian of a few massive equations.
Per the documentation I have set up the problem as follows:
% Define symbols to be used throughout.
syms lat lon u v w h phi theta psi p q r delta1 delta2 delta3 delta4 Hg GM J2 Re e Ixx Ixz Iyy Izz omega m real
% Rotational rate in body axis.
omegabnb = [p;q;r];
% Inertia matrix.
I = [Ixx,0,Ixz;0,Iyy,0;-Ixz,0,Izz];
% Rotational rate of Earth.
omegaeie = [0;0;omega];
% Equatorial radius of curvature.
RE = Re/sqrt(1-e^2*sin(lat)^2);
% Polar radius of curvature.
RN = (Re*(1-e^2))/((1-e^2*sin(lat)^2)^(3/2));
% Euler orientations.
roll = [1,0,0;0,cos(phi),-sin(phi);0,sin(phi),cos(phi)];
pitch = [cos(theta),0,sin(theta);0,1,0;-sin(theta),0,cos(theta)];
yaw = [cos(psi),-sin(psi),0;sin(psi),cos(psi),0;0,0,1];
% Euler Matrix.
Cbn = simplify(yaw*pitch*roll);
Cnb = Cbn';
% LLA Matrix.
[~,Cne] = nRe(lat,lon);
Cne = simplify(Cne);
Cen = Cne';
% Earth rotation orientation.
[~,Cei] = eRi(Hg);
Cei = simplify(Cei);
Cie = Cei';
% body velocity
vbeb = [u;v;w];
vneb = Cbn'*vbeb;
vneb = simplify(vneb);
% Latitude and longitudinal rates.
latdot = simplify(vneb(1)/(RN+h));
longdot = simplify((vneb(2)/(RE+h))/cos(lat));
omeganen = [longdot*cos(lat);-latdot;-longdot*sin(lat)];
% Forces (outputs six equations determining the various forces and moments, equations in terms of most of the state)
% VERY long equations (126 coefficients).
[Fb,Mb] = MomentFunc;
% Radii calculations.
Rc = [1/(RN+h),0,0;
0,1/((RN+h)*cos(lat)),0 ;
0,0,-1];
% Position calculation.
reeb = [((Re/sqrt(1-e^2*(sin(lat)^2)))+h)*cos(lat)*cos(lon);
((Re/sqrt(1-e^2*(sin(lat)^2)))+h)*cos(lat)*sin(lon);
((Re*(1-e^2)/sqrt(1-e^2*(sin(lat)^2)))+h)*sin(lat)];
riib = Cie*reeb;
Riib = sqrt(riib(1)^2+riib(2)^2+riib(3)^2);
tmpx = (1-5*(riib(3)/Riib)^2)*riib(1);
tmpy = (1-5*(riib(3)/Riib)^2)*riib(2);
tmpz = (3-5*(riib(3)/Riib)^2)*riib(3);
tmp = [tmpx;tmpy;tmpz];
tmp = riib + (3/2)*J2*(Re^2)*tmp/(Riib^2);
% Acceleration due to gravity.
gammaiib = -GM/(Riib^3)*(tmp);
% Latitude, Longitude and Altitude rate.
gammadot = Rc*vneb;
% Body velocity rate.
vbebdot = Cen*Cnb*(Fb/m) + Cei*gammaiib - Cbn*Cne*SkewMat(Cen*omeganen+Cne*Cnb*omegabnb)*Cen*Cnb*vbeb;
% Angular rotational rate.
omegabnbdot = I\(Mb-SkewMat(Cbn*Cne*omegaeie+Cbn*omeganen+omegabnb)*I*(Cbn*Cne*omegaeie+Cbn*omeganen+omegabnb))-SkewMat(-omegabnb)*Cbn*omeganen;
% Roll rate.
phidot = p*q*sin(phi)*tan(theta)+r*cos(phi)*tan(theta);
% Pitch rate.
thetadot = q*cos(phi)-r*sin(phi);
% Yaw rate.
psidot = q*sin(phi)*sec(theta)+r*cos(phi)*sec(theta);
% Control deflection rate (DUE TO STATE, NOT DUE TO CONTROL).
deltadot = 0;
% Jacobian for linearization.
A = jacobian([gammadot;vbebdot;omegabnbdot;phidot;thetadot;psidot;deltadot;deltadot;deltadot;deltadot],[lat,lon,h,u,v,w,p,q,r,phi,theta,psi,delta1,delta2,delta3,delta4]);
This progam runs for several days and then ends with an out of memory error. I wish I'd stored the error, but I will have it once it finishes running again.
For reference I have 64GB ram.
Unclear above is the fact that vbebdot and omegabnbdot are HUGE equations (spanning several pages).

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!