Solving nonLinear equation and Linearized Eqaution using ODE45

14 vues (au cours des 30 derniers jours)
Alex
Alex le 9 Déc 2019
Commenté : darova le 10 Déc 2019
I've been doing a comparison of the plots for roll, yaw, and pitch of a satellite in GG orbit and I've been unable to get my nonLinear function to relatively match my Linearized equation. I've tried my variant of the equations and the provided variants while only getting the linearized equation to come out correctly. As the linearized equation is just a simplified version of the nonlinear equation, and I've tried known good equations without success, I suspect I am doing something wrong in how I set up the code.
Here is the base equation's that are being worked off off of. The first longer one is set equal to the smaller second one to make the end result that is used in the code for the nonlinear function. This is the same equation set that when linearized, works just fine in ODE45. Since the first one is so long, this is why I broke it down into a smaller equation by using alphabetical variable for the solver.
Equations_long.png
Equation_short.png
To start off, the three equations are dependent on each other a little bit so I simplified the equations to alphabetical variables, and then had MATLAB solve them for me. I then just add in what the variables are equal to once in the the ODE45 function. I must be missing something here. Below is the code for my derivation of the equations, plotting results are similar to the know correct ones.
% Simplified nonLinear equations
eq1 = ((A-B+C-D*x6d)/(-F))-x2d;
eq2 = ((G+H*I+L*x6d)/(-K))-x4d;
eq3 = ((M-N*O-P*x4d)/(-R))-x6d;
% Solve 2nd order ODE's to get in terms of x2d, x4d, and x6d
[sol1, sol2, sol3] = solve([eq1, eq2, eq3],[x2d,x4d,x6d]);
% Initial conditions of 1 degree in yaw, roll, pitch
X0 = [1 0 1 0 1 0];
% 1.6 times the rotation period of satellite in seconds
t = 0:60:5676.98*1.6;
[time_lin, solution_lin] = ode45(@(t,x) Lin_func(t,x), t, X0);
[time_nonlin, solution_nonlin] = ode45(@(t,x) nonLin_func(t,x), t, X0);
% X has [X(1), X(2), X(3), X(4), X(5), X(6)]
function d_dt = Lin_func(t,X)
% Constants
Ix = 6;
Iy = 8;
Iz = 4;
mu = (3.986*10^14);
Rad = 6878140;
w = sqrt(mu/Rad^3);
% Generate ODE's
d_dt = X*NaN;
d_dt(1) = X(2);
d_dt(2) = (w*(Ix+Iz-Iy)*X(6)-4*(w^2)*(Iy-Iz)*X(1))/Ix;
d_dt(3) = X(4);
d_dt(4) = 3*(w^2)*(Iz-Ix)*X(3)/Iy;
d_dt(5) = X(6);
d_dt(6) = ((w^2)*(Ix-Iy)*X(5)-w*(Iz+Ix-Iy)*X(2))/Iz;
end
function d_dt = nonLin_func(t,X)
% Constants
Ix = 6;
Iy = 8;
Iz = 4;
mu = (3.986*10^14);
Rad = 6878140;
w = sqrt(mu/Rad^3);
% Equation chunks
A = (Iy-Iz)*(w*cosd(X(5))*sind(X(1))-w*cosd(X(1))*sind(X(3))*sind(X(5))-sind(X(1))*X(4)+cosd(X(3))*cosd(X(1))*X(6))*(w*cosd(X(1))*cosd(X(5))+w*sind(X(3))*sind(X(1))*sind(X(5))-cosd(X(1))*X(4)-cosd(X(3))*sind(X(1))*X(6));
B = Ix*w*cosd(X(3))*cosd(X(5))*X(6);
C = Ix*X(4)*(w*sind(X(3))*sind(X(5))-cosd(X(3))*X(6));
D = Ix*sind(X(3));
E = (3*mu/(2*Rad^3))*(sind(2*X(1))*(cosd(X(3))^2)*(Iy-Iz));
F = Ix;
G = (Iz-Ix)*(w*cosd(X(5))*sind(X(1))-w*cosd(X(1))*sind(X(3))*sind(X(5))-sind(X(1))*X(4)+cosd(X(3))*cosd(X(1))*X(6))*(w*cosd(X(3))*sind(X(5))-X(2)+sind(X(3))*X(6));
H = Iy;
I = (-w*cosd(X(5))*sind(X(3))*sind(X(1))*X(6)+w*cosd(X(1))*sind(X(5))*X(6)+X(2)*(w*cosd(X(5))*sind(X(1))-w*cosd(X(1))*sind(X(3))*sind(X(5))+cosd(X(3))*cosd(X(1))*X(6))-sind(X(1))*X(4)*(w*cosd(X(3))*sind(X(5))+X(2)+sind(X(3))*X(6)));
J = (3*mu/(2*Rad^3))*(sind(2*X(3))*cosd(X(1))*(Ix-Iz));
K = cosd(X(1))*Iy;
L = Iy*cosd(X(3))*sind(X(1));
M = (Iy-Ix)*(-w*cosd(X(3))*sind(X(5))+X(2)-sind(X(3))*X(6))*(-w*(cosd(X(1))*cosd(X(5))+sind(X(3))*sind(X(1))*sind(X(5)))+cosd(X(1))*X(4)+cosd(X(3))*sind(X(1))*X(6));
N = Iz;
O = (w*cosd(X(1))*cosd(X(5))*sind(X(3))*X(6)+w*sind(X(1))*sind(X(5))*X(6)+cosd(X(1))*X(4)*(w*cosd(X(3))*sind(X(5))+X(2)+sind(X(3))*X(6))-X(2)*(w*cosd(X(1))*cosd(X(5))+w*sind(X(3))*sind(X(1))*sind(X(5))-cosd(X(3))*sind(X(1))*X(6)));
P = sind(X(1))*Iz;
Q = (3*mu/(2*Rad^3))*(sind(2*X(3))*sind(X(1))*(Iy-Ix));
R = cosd(X(3))*cosd(X(1))*Iz;
% Generate ODE's
d_dt = X*NaN;
d_dt(1) = X(2);
d_dt(2) = -(D*G*P + D*K*M + A*L*P + A*K*R - B*L*P - D*J*P - B*K*R + C*L*P + C*K*R - D*K*Q - E*L*P - E*K*R + D*H*I*P - D*K*N*O)/(F*(L*P + K*R));
d_dt(3) = X(4);
d_dt(4) = -(G*R - L*M - J*R + L*Q + H*I*R + L*N*O)/(L*P + K*R);
d_dt(5) = X(6);
d_dt(6) = -(G*P + K*M - J*P - K*Q + H*I*P - K*N*O)/(L*P + K*R);
end
  2 commentaires
darova
darova le 9 Déc 2019
How to check? Where are original equations/formujlas?
Alex
Alex le 9 Déc 2019
I apologize for leaving out the original equations that the code is based off of. I have updated my posting to include these.

Connectez-vous pour commenter.

Réponse acceptée

darova
darova le 9 Déc 2019
YOu have to express , and . Solve the system
image.png
function du = nonlin(t,y)
fi = y(1);
psi = y(2);
theta = y(3);
% shortcuts
cf = cosd(fi);
cp = cosd(psi);
ct = cosd(theta);
sf = sind(fi);
sp = sind(psi);
st = sind(theta);
% coefficients of fi'', psi'' and theta''
A = [1 -st 0
0 ct*sf cf
0 ct*cf sf];
B(1,1) = long expression1;
B(2,1) = long expression2;
B(3,1) = long expression3;
X = A\B;
du(1:3,1) = y(4:6);
du(4:6,1) = X;
end
  2 commentaires
Alex
Alex le 10 Déc 2019
Do I input the three expressions as is, or do I input them in solved for fi'', psi'' and theta''?
darova
darova le 10 Déc 2019
I highlighted long expressions. Just change the sign
image.png

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Agriculture dans Help Center et File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by