ode45 problems -- not getting expected solution

25 次查看（过去 30 天）
Rick Sellers
Rick Sellers 2020-7-22
编辑： David Goodmanson 2020-7-24
Hello,
I have a set of nonlinear satellite attitude equations to solve. I've implemented this as I've done in the past (at least I think I've done everything the same), and the solutions I get are linear. There should be an oscillation in each of the angles phi, tht, psi over time. My suspicion is that I didn't define the vector correctly being sent to ode45, my ode function doesn't update the contents of this vector correctly or there's some sort of scaling issue. I've asked other classmates already, the TA for my class, and stared at it for hours without seeing the error, so I'm turning to the experts for assistance.
Thank you,
Rick
% Derive Equations of Motion
clc; clear;
% PRY = [phi,tht,psi,dphi,dtht,dpsi]';
% Create sym variables
syms Sphi(t) Stht(t) Spsi(t) dSphi(t) dStht(t) dSpsi(t) ddSphi(t) ddStht(t)...
ddSpsi(t)
assume([Sphi(t) Stht(t) Spsi(t) dSphi(t) dStht(t) dSpsi(t) ddSphi(t) ...
ddStht(t) ddSpsi(t)],'real')
syms SOmega SR SG SmE SIxx SIyy SIzz real positive
% Form rotation matrices and cosine direction matrices
Rxphi = [1,0,0;0,cos(Sphi),sin(Sphi);0,-sin(Sphi),cos(Sphi)];
Rytht(t) = [cos(Stht),0,-sin(Stht);0,1,0;sin(Stht),0,cos(Stht)];
Rzpsi(t) = [cos(Spsi),sin(Spsi),0;-sin(Spsi),cos(Spsi),0;0,0,1];
BCA = Rxphi*Rytht*Rzpsi;
ACB = transpose(BCA);
% Compute O_omega_A in B-frame
OwA_B = BCA*[0;0;SOmega];
% Compute A_omega_B in B-frame
AwBcross = BCA*diff(ACB,t);
AwB_Bx = [0,0,1]*AwBcross*[0;1;0];
AwB_By = [1,0,0]*AwBcross*[0;0;1];
AwB_Bz = [0,1,0]*AwBcross*[1;0;0];
AwB_B = [ AwB_Bx;AwB_By;AwB_Bz ];
% Compute O_omega_B in B-frame
OwB_B = OwA_B + AwB_B;
% Decompose O_omega_B
OwB = eye(3)*OwB_B(t);
OwB = [OwB(1,:);OwB(2,:);OwB(3,:)];
% Moment of Inertia Matrix
I = [SIxx,0,0;0,SIyy,0;0,0,SIzz];
% Angular Momentum
h = I*OwB;
hdot = diff(h,t);
% Torques
K = 3*SG*SmE/SR^5;
Gx = K*(sin(2*Sphi)*(cos(Stht))^2*(SIzz-SIyy));
Gy = K*(sin(2*Stht)*cos(Sphi)*(SIzz-SIxx));
Gz = K*(sin(2*Stht)*sin(Sphi)*(SIxx-SIyy));
Tau = hdot + cross(OwB,h);
TauX = Tau(1,:);
TauY = Tau(2,:);
TauZ = Tau(3,:);
% Substitute
TauX = subs(TauX,[diff(Sphi(t), t, t),diff(Stht(t), t, t),diff(Spsi(t), t, t)],[ddSphi,ddStht,ddSpsi]);
TauY = subs(TauY,[diff(Sphi(t), t, t),diff(Stht(t), t, t),diff(Spsi(t), t, t)],[ddSphi,ddStht,ddSpsi]);
TauZ = subs(TauZ,[diff(Sphi(t), t, t),diff(Stht(t), t, t),diff(Spsi(t), t, t)],[ddSphi,ddStht,ddSpsi]);
TauX = subs(TauX,[diff(Sphi(t), t),diff(Stht(t), t),diff(Spsi(t), t)],[dSphi,dStht,dSpsi]);
TauY = subs(TauY,[diff(Sphi(t), t),diff(Stht(t), t),diff(Spsi(t), t)],[dSphi,dStht,dSpsi]);
TauZ = subs(TauZ,[diff(Sphi(t), t),diff(Stht(t), t),diff(Spsi(t), t)],[dSphi,dStht,dSpsi]);
TauX = subs(TauX,[Sphi(t),Stht(t),Spsi(t)],[Sphi,Stht,Spsi]);
TauY = subs(TauY,[Sphi(t),Stht(t),Spsi(t)],[Sphi,Stht,Spsi]);
TauZ = subs(TauZ,[Sphi(t),Stht(t),Spsi(t)],[Sphi,Stht,Spsi]);
% Tx = TauX == Gx;
% Ty = TauY == Gy;
% Tz = TauZ == Gz;
Tx = TauX == 0;
Ty = TauY == 0;
Tz = TauZ == 0;
isolate(Tx,ddSphi)
isolate(Ty,ddStht)
isolate(Tz,ddSpsi)
% Main
% Sets intitial values, creates PRY vector of unknowns and paramters; calls
% other functions and solver; sorts data and plots results
clc; clear; close all
% Set initial values
PRY = zeros(6,1);
phi = 0; tht = 0; psi = 0;
dphi = 0; dtht = 0; dpsi = 0;
RE = 6378.137e3 ; % m - radius of earth
RS = 500e3; % m - orbital altitude
R = (RE+RS); % m
mE = 5.972e24; % kg
G = 6.6743e-11; % m^3/kg s^2
Omega = sqrt(G*mE/R^3); % rad/s
Ixx = 6; Iyy = 8; Izz = 4; % kg m^2
mu = G*mE;
per = 2*pi*(sqrt(R^3/mu));
tend = 1.6*per;
tspan = [0 tend];
ic = [1;5;10;15;20;25];
for k= 1:6
phi =ic(k);
tht=ic(k);
psi=ic(k);
PRY = [phi,tht,psi,dphi,dtht,dpsi]';
[t,prydot] = ode45(@(t,pry) proj_ODE(t,PRY,Omega,R,mE,G,Ixx,Iyy,Izz),tspan,PRY);
% [t,prydot] = ode45(@(t,pry) proj_ODE_L(t,PRY,Omega,R,mE,G,Ixx,Iyy,Izz),tspan,PRY);
% Plotting results; pitch, roll, yaw for each initial condition
xticks = linspace(0,1.6,9);
figure(k)
% prydot = rad2deg(prydot);
subplot(3,1,1)
plot(t,prydot(:,1))
grid on
title('Pitch, \phi');
xlabel('Time, [s]'); ylabel('');
subplot(3,1,2)
plot(t,prydot(:,2));
grid on
title('Roll, \theta');
xlabel('Time, [s]');
subplot(3,1,3)
plot(t,prydot(:,3));
grid on
title('Yaw, \psi');
xlabel('Time, [s]'); ylabel('Attitude, [rad]');
end
% pry_ODE.m
% myODE function for solver
% PRY = [phi,tht,psi,dphi,dtht,dpsi];
function prydot = proj_ODE(t,pry,Omega,R,mE,G,Ixx,Iyy,Izz)
prydot = zeros(size(pry));
phi = pry(1);
tht = pry(2);
psi = pry(3);
dphi = pry(4);
dtht = pry(5);
dpsi = pry(6);
K = 3*G*mE/R^5;
w0 = Omega;
Gx = K*(sin(2*phi)*(cos(tht))^2*(Izz-Iyy));
Gy = K*(sin(2*tht)*cos(phi)*(Izz-Ixx));
Gz = K*(sin(2*tht)*sin(phi)*(Ixx-Iyy));
prydot(1) = dphi;
prydot(2) = dtht;
prydot(3) = dpsi;
prydot(4) = (Ixx*((sin(phi)*sin(psi) ...
+ cos(phi)*cos(psi)*sin(tht))*(cos(phi)*cos(psi)*prydot(6) ...
- cos(phi)*sin(psi)*dphi^2 - cos(phi)*sin(psi)*dpsi^2 ...
+ cos(psi)*sin(phi)*sin(tht)*dphi^2 + cos(psi)*sin(phi)*sin(tht)*dpsi^2 ...
+ cos(psi)*sin(phi)*sin(tht)*dtht^2 - 2*cos(psi)*sin(phi)*dphi*dpsi ...
- cos(psi)*cos(tht)*sin(phi)*prydot(5) + sin(phi)*sin(psi)*sin(tht)*prydot(6) ...
- 2*cos(phi)*cos(psi)*cos(tht)*dphi*dtht ...
+ 2*cos(phi)*sin(psi)*sin(tht)*dphi*dpsi ...
+ 2*cos(tht)*sin(phi)*sin(psi)*dpsi*dtht) - (cos(psi)*sin(phi) ...
- cos(phi)*sin(psi)*sin(tht))*(cos(phi)*sin(psi)*prydot(6) ...
+ cos(phi)*cos(psi)*dphi^2 + cos(phi)*cos(psi)*dpsi^2 ...
+ sin(phi)*sin(psi)*sin(tht)*dphi^2 + sin(phi)*sin(psi)*sin(tht)*dpsi^2 ...
+ sin(phi)*sin(psi)*sin(tht)*dtht^2 - 2*sin(phi)*sin(psi)*dphi*dpsi ...
- cos(psi)*sin(phi)*sin(tht)*prydot(6) - cos(tht)*sin(phi)*sin(psi)*prydot(5) ...
- 2*cos(phi)*cos(psi)*sin(tht)*dphi*dpsi - 2*cos(phi)*cos(tht)*sin(psi)*dphi*dtht ...
- 2*cos(psi)*cos(tht)*sin(phi)*dpsi*dtht) - (cos(phi)*sin(psi)*dphi ...
+ cos(psi)*sin(phi)*dpsi + cos(phi)*cos(psi)*cos(tht)*dtht ...
- cos(psi)*sin(phi)*sin(tht)*dphi ...
- cos(phi)*sin(psi)*sin(tht)*dpsi)*(sin(phi)*sin(psi)*dphi ...
- cos(phi)*cos(psi)*dpsi + cos(phi)*cos(psi)*sin(tht)*dphi ...
+ cos(psi)*cos(tht)*sin(phi)*dtht - sin(phi)*sin(psi)*sin(tht)*dpsi) ...
- (sin(phi)*sin(psi)*dpsi - cos(phi)*cos(psi)*dphi ...
+ cos(phi)*cos(psi)*sin(tht)*dpsi + cos(phi)*cos(tht)*sin(psi)*dtht ...
- sin(phi)*sin(psi)*sin(tht)*dphi)*(cos(phi)*sin(psi)*sin(tht)*dphi ...
- cos(phi)*sin(psi)*dpsi - cos(psi)*sin(phi)*dphi ...
+ cos(psi)*sin(phi)*sin(tht)*dpsi + cos(tht)*sin(phi)*sin(psi)*dtht) ...
+ cos(phi)*cos(tht)*(sin(phi)*sin(tht)*prydot(5) ...
+ cos(tht)*sin(phi)*dphi^2 + cos(tht)*sin(phi)*dtht^2 ...
+ 2*cos(phi)*sin(tht)*dphi*dtht) + w0*cos(tht)*dtht ...
+ cos(tht)*sin(phi)*dphi*(cos(phi)*cos(tht)*dphi - sin(phi)*sin(tht)*dtht) ...
+ cos(phi)*sin(tht)*dtht*(cos(phi)*cos(tht)*dphi - sin(phi)*sin(tht)*dtht)) ...
+ Iyy*(sin(tht)*(cos(tht)*sin(phi)*dphi + cos(phi)*sin(tht)*dtht) ...
+ w0*cos(tht)*sin(phi) + cos(psi)*cos(tht)*(cos(phi)*sin(psi)*dphi ...
+ cos(psi)*sin(phi)*dpsi + cos(phi)*cos(psi)*cos(tht)*dtht ...
- cos(psi)*sin(phi)*sin(tht)*dphi - cos(phi)*sin(psi)*sin(tht)*dpsi) ...
+ cos(tht)*sin(psi)*(sin(phi)*sin(psi)*dpsi - cos(phi)*cos(psi)*dphi ...
+ cos(phi)*cos(psi)*sin(tht)*dpsi + cos(phi)*cos(tht)*sin(psi)*dtht ...
- sin(phi)*sin(psi)*sin(tht)*dphi))*((cos(phi)*cos(psi) ...
+ sin(phi)*sin(psi)*sin(tht))*(cos(psi)*cos(tht)*dpsi ...
- sin(psi)*sin(tht)*dtht) + (cos(phi)*sin(psi) ...
- cos(psi)*sin(phi)*sin(tht))*(cos(tht)*sin(psi)*dpsi ...
+ cos(psi)*sin(tht)*dtht) + w0*cos(phi)*cos(tht) - cos(tht)^2*sin(phi)*dtht) ...
- Izz*(sin(tht)*(cos(tht)*sin(phi)*dphi + cos(phi)*sin(tht)*dtht) ...
+ w0*cos(tht)*sin(phi) + cos(psi)*cos(tht)*(cos(phi)*sin(psi)*dphi ...
+ cos(psi)*sin(phi)*dpsi + cos(phi)*cos(psi)*cos(tht)*dtht ...
- cos(psi)*sin(phi)*sin(tht)*dphi - cos(phi)*sin(psi)*sin(tht)*dpsi) ...
+ cos(tht)*sin(psi)*(sin(phi)*sin(psi)*dpsi - cos(phi)*cos(psi)*dphi ...
+ cos(phi)*cos(psi)*sin(tht)*dpsi + cos(phi)*cos(tht)*sin(psi)*dtht ...
- sin(phi)*sin(psi)*sin(tht)*dphi))*((cos(phi)*cos(psi) ...
+ sin(phi)*sin(psi)*sin(tht))*(cos(psi)*cos(tht)*dpsi ...
- sin(psi)*sin(tht)*dtht) + (cos(phi)*sin(psi) ...
- cos(psi)*sin(phi)*sin(tht))*(cos(tht)*sin(psi)*dpsi ...
+ cos(psi)*sin(tht)*dtht) + w0*cos(phi)*cos(tht) ...
- cos(tht)^2*sin(phi)*dtht))/(Ixx*(cos(phi)^2*cos(tht)^2 + (sin(phi)*sin(psi) ...
+ cos(phi)*cos(psi)*sin(tht))^2 + (cos(psi)*sin(phi) ...
- cos(phi)*sin(psi)*sin(tht))^2));
prydot(5) = (Iyy*(cos(tht)*sin(psi)*(cos(phi)*cos(psi)*prydot(4) ...
- sin(phi)*sin(psi)*prydot(6) - cos(psi)*sin(phi)*dphi^2 ...
- cos(psi)*sin(phi)*dpsi^2 + cos(phi)*sin(psi)*sin(tht)*dphi^2 ...
+ cos(phi)*sin(psi)*sin(tht)*dpsi^2 + cos(phi)*sin(psi)*sin(tht)*dtht^2 ...
- 2*cos(phi)*sin(psi)*dphi*dpsi - cos(phi)*cos(psi)*sin(tht)*prydot(6) ...
+ sin(phi)*sin(psi)*sin(tht)*prydot(4) ...
- 2*cos(phi)*cos(psi)*cos(tht)*dpsi*dtht ...
+ 2*cos(psi)*sin(phi)*sin(tht)*dphi*dpsi ...
+ 2*cos(tht)*sin(phi)*sin(psi)*dphi*dtht) ...
- sin(tht)*(cos(tht)*sin(phi)*prydot(4) + cos(phi)*cos(tht)*dphi^2 ...
+ cos(phi)*cos(tht)*dtht^2 - 2*sin(phi)*sin(tht)*dphi*dtht) ...
+ cos(psi)*cos(tht)*(sin(phi)*sin(psi)*dphi^2 + sin(phi)*sin(psi)*dpsi^2 ...
- cos(phi)*sin(psi)*prydot(4) - cos(psi)*sin(phi)*prydot(6) ...
+ cos(phi)*cos(psi)*sin(tht)*dphi^2 + cos(phi)*cos(psi)*sin(tht)*dpsi^2 ...
+ cos(phi)*cos(psi)*sin(tht)*dtht^2 - 2*cos(phi)*cos(psi)*dphi*dpsi ...
+ cos(psi)*sin(phi)*sin(tht)*prydot(4) ...
+ cos(phi)*sin(psi)*sin(tht)*prydot(6) ...
+ 2*cos(psi)*cos(tht)*sin(phi)*dphi*dtht ...
+ 2*cos(phi)*cos(tht)*sin(psi)*dpsi*dtht ...
- 2*sin(phi)*sin(psi)*sin(tht)*dphi*dpsi) ...
- cos(tht)*dtht*(cos(tht)*sin(phi)*dphi + cos(phi)*sin(tht)*dtht) ...
+ w0*sin(phi)*sin(tht)*dtht - cos(psi)*cos(tht)*dpsi*(sin(phi)*sin(psi)*dpsi ...
- cos(phi)*cos(psi)*dphi + cos(phi)*cos(psi)*sin(tht)*dpsi ...
+ cos(phi)*cos(tht)*sin(psi)*dtht - sin(phi)*sin(psi)*sin(tht)*dphi) ...
+ cos(tht)*sin(psi)*dpsi*(cos(phi)*sin(psi)*dphi + cos(psi)*sin(phi)*dpsi ...
+ cos(phi)*cos(psi)*cos(tht)*dtht - cos(psi)*sin(phi)*sin(tht)*dphi ...
- cos(phi)*sin(psi)*sin(tht)*dpsi) ...
+ cos(psi)*sin(tht)*dtht*(cos(phi)*sin(psi)*dphi ...
+ cos(psi)*sin(phi)*dpsi + cos(phi)*cos(psi)*cos(tht)*dtht ...
- cos(psi)*sin(phi)*sin(tht)*dphi - cos(phi)*sin(psi)*sin(tht)*dpsi) ...
+ sin(psi)*sin(tht)*dtht*(sin(phi)*sin(psi)*dpsi - cos(phi)*cos(psi)*dphi ...
+ cos(phi)*cos(psi)*sin(tht)*dpsi + cos(phi)*cos(tht)*sin(psi)*dtht ...
- sin(phi)*sin(psi)*sin(tht)*dphi) - w0*cos(phi)*cos(tht)*dphi) ...
- Ixx*((cos(phi)*cos(psi) + sin(phi)*sin(psi)*sin(tht))*(cos(psi)*cos(tht)*dpsi ...
- sin(psi)*sin(tht)*dtht) + (cos(phi)*sin(psi) ...
- cos(psi)*sin(phi)*sin(tht))*(cos(tht)*sin(psi)*dpsi ...
+ cos(psi)*sin(tht)*dtht) + w0*cos(phi)*cos(tht) ...
- cos(tht)^2*sin(phi)*dtht)*((sin(phi)*sin(psi) ...
+ cos(phi)*cos(psi)*sin(tht))*(sin(phi)*sin(psi)*dphi ...
- cos(phi)*cos(psi)*dpsi + cos(phi)*cos(psi)*sin(tht)*dphi ...
+ cos(psi)*cos(tht)*sin(phi)*dtht - sin(phi)*sin(psi)*sin(tht)*dpsi) ...
- (cos(psi)*sin(phi) ...
- cos(phi)*sin(psi)*sin(tht))*(cos(phi)*sin(psi)*sin(tht)*dphi ...
- cos(phi)*sin(psi)*dpsi - cos(psi)*sin(phi)*dphi ...
+ cos(psi)*sin(phi)*sin(tht)*dpsi + cos(tht)*sin(phi)*sin(psi)*dtht) ...
- w0*sin(tht) + cos(phi)*cos(tht)*(cos(phi)*cos(tht)*dphi ...
- sin(phi)*sin(tht)*dtht)) + Izz*((cos(phi)*cos(psi) ...
+ sin(phi)*sin(psi)*sin(tht))*(cos(psi)*cos(tht)*dpsi ...
- sin(psi)*sin(tht)*dtht) + (cos(phi)*sin(psi) ...
- cos(psi)*sin(phi)*sin(tht))*(cos(tht)*sin(psi)*dpsi ...
+ cos(psi)*sin(tht)*dtht) + w0*cos(phi)*cos(tht) ...
- cos(tht)^2*sin(phi)*dtht)*((sin(phi)*sin(psi) ...
+ cos(phi)*cos(psi)*sin(tht))*(sin(phi)*sin(psi)*dphi ...
- cos(phi)*cos(psi)*dpsi + cos(phi)*cos(psi)*sin(tht)*dphi ...
+ cos(psi)*cos(tht)*sin(phi)*dtht - sin(phi)*sin(psi)*sin(tht)*dpsi) ...
- (cos(psi)*sin(phi) ...
- cos(phi)*sin(psi)*sin(tht))*(cos(phi)*sin(psi)*sin(tht)*dphi ...
- cos(phi)*sin(psi)*dpsi - cos(psi)*sin(phi)*dphi ...
+ cos(psi)*sin(phi)*sin(tht)*dpsi + cos(tht)*sin(phi)*sin(psi)*dtht) ...
- w0*sin(tht) + cos(phi)*cos(tht)*(cos(phi)*cos(tht)*dphi ...
- sin(phi)*sin(tht)*dtht)))/(Iyy*(cos(phi)*cos(psi)^2*cos(tht)^2 ...
+ cos(phi)*cos(tht)^2*sin(psi)^2 + cos(phi)*sin(tht)^2));
prydot(6) = (Izz*((cos(phi)*cos(psi) ...
+ sin(phi)*sin(psi)*sin(tht))*(sin(psi)*sin(tht)*prydot(5) ...
+ cos(tht)*sin(psi)*dpsi^2 + cos(tht)*sin(psi)*dtht^2 ...
+ 2*cos(psi)*sin(tht)*dpsi*dtht) + (cos(tht)*sin(psi)*dpsi ...
+ cos(psi)*sin(tht)*dtht)*(sin(phi)*sin(psi)*dphi ...
- cos(phi)*cos(psi)*dpsi + cos(phi)*cos(psi)*sin(tht)*dphi ...
+ cos(psi)*cos(tht)*sin(phi)*dtht - sin(phi)*sin(psi)*sin(tht)*dpsi) ...
- (cos(psi)*cos(tht)*dpsi ...
- sin(psi)*sin(tht)*dtht)*(cos(phi)*sin(psi)*sin(tht)*dphi ...
- cos(phi)*sin(psi)*dpsi - cos(psi)*sin(phi)*dphi ...
+ cos(psi)*sin(phi)*sin(tht)*dpsi + cos(tht)*sin(phi)*sin(psi)*dtht) ...
- (cos(phi)*sin(psi) ...
- cos(psi)*sin(phi)*sin(tht))*(cos(psi)*sin(tht)*prydot(5) ...
+ cos(psi)*cos(tht)*dpsi^2 + cos(psi)*cos(tht)*dtht^2 ...
- 2*sin(psi)*sin(tht)*dpsi*dtht) + cos(tht)^2*sin(phi)*prydot(5) ...
- 2*cos(tht)*sin(phi)*sin(tht)*dtht^2 + cos(phi)*cos(tht)^2*dphi*dtht ...
+ w0*cos(tht)*sin(phi)*dphi + w0*cos(phi)*sin(tht)*dtht) ...
+ Ixx*(sin(tht)*(cos(tht)*sin(phi)*dphi + cos(phi)*sin(tht)*dtht) ...
+ w0*cos(tht)*sin(phi) + cos(psi)*cos(tht)*(cos(phi)*sin(psi)*dphi ...
+ cos(psi)*sin(phi)*dpsi + cos(phi)*cos(psi)*cos(tht)*dtht ...
- cos(psi)*sin(phi)*sin(tht)*dphi - cos(phi)*sin(psi)*sin(tht)*dpsi) ...
+ cos(tht)*sin(psi)*(sin(phi)*sin(psi)*dpsi - cos(phi)*cos(psi)*dphi ...
+ cos(phi)*cos(psi)*sin(tht)*dpsi + cos(phi)*cos(tht)*sin(psi)*dtht ...
- sin(phi)*sin(psi)*sin(tht)*dphi))*((sin(phi)*sin(psi) ...
+ cos(phi)*cos(psi)*sin(tht))*(sin(phi)*sin(psi)*dphi ...
- cos(phi)*cos(psi)*dpsi + cos(phi)*cos(psi)*sin(tht)*dphi ...
+ cos(psi)*cos(tht)*sin(phi)*dtht - sin(phi)*sin(psi)*sin(tht)*dpsi) ...
- (cos(psi)*sin(phi) ...
- cos(phi)*sin(psi)*sin(tht))*(cos(phi)*sin(psi)*sin(tht)*dphi ...
- cos(phi)*sin(psi)*dpsi - cos(psi)*sin(phi)*dphi ...
+ cos(psi)*sin(phi)*sin(tht)*dpsi + cos(tht)*sin(phi)*sin(psi)*dtht) ...
- w0*sin(tht) + cos(phi)*cos(tht)*(cos(phi)*cos(tht)*dphi ...
- sin(phi)*sin(tht)*dtht)) - Iyy*(sin(tht)*(cos(tht)*sin(phi)*dphi ...
+ cos(phi)*sin(tht)*dtht) + w0*cos(tht)*sin(phi) ...
+ cos(psi)*cos(tht)*(cos(phi)*sin(psi)*dphi + cos(psi)*sin(phi)*dpsi ...
+ cos(phi)*cos(psi)*cos(tht)*dtht - cos(psi)*sin(phi)*sin(tht)*dphi ...
- cos(phi)*sin(psi)*sin(tht)*dpsi) ...
+ cos(tht)*sin(psi)*(sin(phi)*sin(psi)*dpsi - cos(phi)*cos(psi)*dphi ...
+ cos(phi)*cos(psi)*sin(tht)*dpsi + cos(phi)*cos(tht)*sin(psi)*dtht ...
- sin(phi)*sin(psi)*sin(tht)*dphi))*((sin(phi)*sin(psi) ...
+ cos(phi)*cos(psi)*sin(tht))*(sin(phi)*sin(psi)*dphi ...
- cos(phi)*cos(psi)*dpsi + cos(phi)*cos(psi)*sin(tht)*dphi ...
+ cos(psi)*cos(tht)*sin(phi)*dtht - sin(phi)*sin(psi)*sin(tht)*dpsi) ...
- (cos(psi)*sin(phi) ...
- cos(phi)*sin(psi)*sin(tht))*(cos(phi)*sin(psi)*sin(tht)*dphi ...
- cos(phi)*sin(psi)*dpsi - cos(psi)*sin(phi)*dphi ...
+ cos(psi)*sin(phi)*sin(tht)*dpsi + cos(tht)*sin(phi)*sin(psi)*dtht) ...
- w0*sin(tht) + cos(phi)*cos(tht)*(cos(phi)*cos(tht)*dphi ...
- sin(phi)*sin(tht)*dtht)))/(Izz*(cos(psi)*cos(tht)*(cos(phi)*cos(psi) ...
+ sin(phi)*sin(psi)*sin(tht)) + cos(tht)*sin(psi)*(cos(phi)*sin(psi) ...
- cos(psi)*sin(phi)*sin(tht))));
end
  15 个评论
显示 13更早的评论
Rick Sellers
Rick Sellers 2020-7-24
David,
I ended up getting it working. It doesn't quite match the plots from the paper we were given to re-create, but that paper was basically garbage, and so I'm happy with what I have now. As someone else pointed out, my return vector from the ODE function was zeroized at the beginning of the function. This, along with the problem you noted were the source of my issues.
Another comment had asked about the derivation of these equations, adn that's what the EOM code is about. You're correct, that was my back-end symbolic work to get these equations. They're very simply derived from 1) getting the angular velociy components of the satellite body frame wrt an inertial reference frame; 2) computing the angular momentum as h = [I] * omega, where I is the moment of inertia tensor; 3) the torques on the satellite are then given by {T} = d/dt(h) + omega X h; The LHS of the equations comes from gravity gradient moments (the 3 G m/R^3 terms).
I obviously have a very rudimentary understanding of this, but hope this shed some light on what the problem was.
'
Thank you for the help,
Rick
David Goodmanson
David Goodmanson 2020-7-24
编辑：David Goodmanson 2020-7-24
Hi RIck,
It's good that you got it working. I feel that the code would be better (certainly way more compact) without using syms at all, but if it ain't broke ... One thing I didn't mention is that in your calculation of the prydots, prydot(4) is calculated in terms of the old prydots, but prydot(5) is calculated with the new prydot(4), and prydot(6) is calculated with new prydot(4) and new prydot(5). It's more consistent to calculate the new ones strictly in terms of the old ones. This might make only a small difference numerically, or it could have a measurable effect.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Steven Lord
Steven Lord 2020-7-22
K = 3*G*mE/R^5;
It's been a while since I've studied any physics, but this looks odd to me.
RE = 6378.137e3 ; % m - radius of earth
RS = 500e3; % m - orbital altitude
R = (RE+RS); % m
mE = 5.972e24; % kg
G = 6.6743e-11; % m^3/kg s^2
By dimensional analysis, K has units of or times whatever the unit of 3 is. If that's supposed to be the gravitational force between a 3 kg satellite and the Earth, shouldn't that be over R^2 (which would give Newtons, ) instead of over R^5?
  1 个评论
Rick Sellers
Rick Sellers 2020-7-22
Steven,
Those terms are multiplied by the moments of inertia [kg * m^2], but that still leaves a unit of length^2 missing to make it Newton-m (these are torques). So, I'm thinking it should be over R^3. But, I just made this change, ran it, and got the same results.
Thank you
Rick

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

产品


版本

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by