need help with my MATLAB code. I am trying to numerically integrate two sets of equations. I want the ode45 to output a column vector of size (12,1). I am getting an error that says unable to meet integration tolerances without reducing the step size. Can you fix the code and make sure it works? % Initial conditions mu = 398600; % km^3/s^2 R = 7000; % km I = [500; 600; 800] * 10^-6; % kg*km^2 A = [0;0]; % Initial PRP and velocity vectors lambda = [1/sqrt(3); 1/sqrt(3); 1/sqrt(3)]; theta = pi/60; % rads w = [0; 0; sqrt(mu/R^3)]; t = 0:43200; % sec % Finding EP EP = PRP2EP(lambda, theta); % Using ode45 to integrate KDE and EOM options = odeset('RelTol',1e-10,'AbsTol',1e-10); [t, y] = ode45(@KDE_C, t, [EP; w; I; A], options); function EP = PRP2EP(lambda, theta) % Finding EP from PRP EP1 = lambda(1)*sind(theta/2); EP2 = lambda(2)*sind(theta/2); EP3 = lambda(3)*sind(theta/2); EP4 = cosd(theta/2); EP = [EP1; EP2; EP3; EP4]; end function dCwdt = KDE_C(~,EPwI) EP = EPwI(1:4); w = EPwI(5:7); I = EPwI(8:10); A = EPwI(11:12); dwdt = zeros(3,1); C11 = 1-2*EP(2,1)^2-2*EP(3,1)^2; C12 = 2*(EP(1,1)*EP(2,1)+EP(3,1)*EP(4,1)); C13 = 2*(EP(1,1)*EP(3,1)-EP(2,1)*EP(4,1)); C21 = 2*(EP(1,1)*EP(2,1)-EP(3,1)*EP(4,1)); C22 = 1-2*EP(1,1)^2-2*EP(3,1)^2; C23 = 2*(EP(3,1)*EP(2,1)+EP(1,1)*EP(4,1)); C31 = 2*(EP(1,1)*EP(3,1)+EP(2,1)*EP(4,1)); C32 = 2*(EP(3,1)*EP(2,1)-EP(1,1)*EP(4,1)); C33 = 1-2*EP(1,1)^2-2*EP(2,1)^2; C11_dot = C12*w(3) - C13*w(2); C12_dot = C13*w(1) - C11*w(3); C13_dot = C11*w(2) - C12*w(1); C21_dot = C22*w(3) - C23*w(2); C22_dot = C23*w(1) - C21*w(3); C23_dot = C21*w(2) - C22*w(1); C31_dot = C32*w(3) - C33*w(2); C32_dot = C33*w(1) - C31*w(3); C33_dot = C31*w(2) - C32*w(1); dCdt1 = [C11_dot; C12_dot; C13_dot]; dCdt2 = [C21_dot; C22_dot; C23_dot]; dCdt3 = [C31_dot; C32_dot; C33_dot]; K1 = (I(2) - I(3)) / I(1); K2 = (I(3) - I(1)) / I(2); K3 = (I(1) - I(2)) / I(3); R_mag = 7000; R1 = R_mag*C11; R2 = R_mag*C21; R3 = R_mag*C31; dwdt(1) = K1*w(2)*w(3) - ((3*w(3)^2*K1*R3*R2)/R_mag^2); dwdt(2) = K2*w(1)*w(3) - ((3*w(3)^2*K2*R1*R3)/R_mag^2); dwdt(3) = K3*w(1)*w(2) - ((3*w(3)^2*K3*R2*R1)/R_mag^2); % Combine the time derivatives into a single vector dCwdt = zeros(12,1); dCwdt = [dCdt1; dCdt2; dCdt3; dwdt]; end
I need help with my MATLAB code. I am trying to numerically integrate two sets of equations. I want the ode45 to output a column vector of size (12,1). I am getting an error that says unable to meet integration tolerances without reducing the step size. Can you fix the code and make sure it works?
% Initial conditions
mu = 398600; % km^3/s^2
R = 7000; % km
I = [500; 600; 800] * 10^-6; % kg*km^2
A = [0;0];
% Initial PRP and velocity
lambda = [1/sqrt(3); 1/sqrt(3); 1/sqrt(3)];
theta = pi/60; % rads
w = [0; 0; sqrt(mu/R^3)];
t = 0:43200; % sec
% Finding EP
EP = PRP2EP(lambda, theta);
% Using ode45 to integrate KDE and EOM
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
[t, y] = ode45(@KDE_C, t, [EP; w; I; A], options);
function EP = PRP2EP(lambda, theta)
% Finding EP from PRP
EP1 = lambda(1)*sind(theta/2);
EP2 = lambda(2)*sind(theta/2);
EP3 = lambda(3)*sind(theta/2);
EP4 = cosd(theta/2);
EP = [EP1; EP2; EP3; EP4];
end
function dCwdt = KDE_C(~,EPwI)
EP = EPwI(1:4);
w = EPwI(5:7);
I = EPwI(8:10);
A = EPwI(11:12);
dwdt = zeros(3,1);
C11 = 1-2*EP(2,1)^2-2*EP(3,1)^2;
C12 = 2*(EP(1,1)*EP(2,1)+EP(3,1)*EP(4,1));
C13 = 2*(EP(1,1)*EP(3,1)-EP(2,1)*EP(4,1));
C21 = 2*(EP(1,1)*EP(2,1)-EP(3,1)*EP(4,1));
C22 = 1-2*EP(1,1)^2-2*EP(3,1)^2;
C23 = 2*(EP(3,1)*EP(2,1)+EP(1,1)*EP(4,1));
C31 = 2*(EP(1,1)*EP(3,1)+EP(2,1)*EP(4,1));
C32 = 2*(EP(3,1)*EP(2,1)-EP(1,1)*EP(4,1));
C33 = 1-2*EP(1,1)^2-2*EP(2,1)^2;
C11_dot = C12*w(3) - C13*w(2);
C12_dot = C13*w(1) - C11*w(3);
C13_dot = C11*w(2) - C12*w(1);
C21_dot = C22*w(3) - C23*w(2);
C22_dot = C23*w(1) - C21*w(3);
C23_dot = C21*w(2) - C22*w(1);
C31_dot = C32*w(3) - C33*w(2);
C32_dot = C33*w(1) - C31*w(3);
C33_dot = C31*w(2) - C32*w(1);
dCdt1 = [C11_dot; C12_dot; C13_dot];
dCdt2 = [C21_dot; C22_dot; C23_dot];
dCdt3 = [C31_dot; C32_dot; C33_dot];
K1 = (I(2) - I(3)) / I(1);
K2 = (I(3) - I(1)) / I(2);
K3 = (I(1) - I(2)) / I(3);
R_mag = 7000;
R1 = R_mag*C11;
R2 = R_mag*C21;
R3 = R_mag*C31;
dwdt(1) = K1*w(2)*w(3) - ((3*w(3)^2*K1*R3*R2)/R_mag^2);
dwdt(2) = K2*w(1)*w(3) - ((3*w(3)^2*K2*R1*R3)/R_mag^2);
dwdt(3) = K3*w(1)*w(2) - ((3*w(3)^2*K3*R2*R1)/R_mag^2);
% Combine the time derivatives into a single vector
dCwdt = zeros(12,1);
dCwdt = [dCdt1; dCdt2; dCdt3; dwdt];
end
Step by step
Solved in 3 steps