How to solve 3 equations dependent each other?

Hello. How to code this situation?

In this situation i need to calculate all angular velocities from time 0s to time 5400s with increase by 0.1s and graph them.

I have all the initial values:

wxi = 0.0065
wyi = 0.0066
wzi = 0.0067

And constant values:

NT = 3.6 * 10^-10
Jx = 2.1 * 10^-3
Jy = 2.0 * 10^-3
Jz = 1.9 * 10^-3

Actually i thought it is going to be easy but when i write the code for wx, i realized that wx is increasing with time but wz and wy in the equation are constant values. Which kind of loop should i use? Thanks for the answers!

[Walter Roberson]

They look to be changing to me. At least one of your J variables would need to be 0 or infinite for a w value to be constant.

[hknatas]

I mean when i calculate wx of course wx is changing with time. But when program calculates such as wx4988, it is always using intial values for wy and wzi but it should use wy4987 and wz4987. I need to a loop such that these equations should be evaluated at the same time from time 0s to time 5400s increase with by 0.1s. Sorry for bad explanation.

[Walter Roberson]

wz(i+1) = wz(i) + delta_t / Jz * (wx(i) * wy(i) + NT)*(Jx - Jy);

[hknatas]

Thanks a lot, it really helps. However i guess i should use a few command which i don't know. As i understand it, i should create wx, wy and wz arrays first. My arrays should have 54000 data and all of them should be equal to their initial values at first. Then after equations are executed one time, second data in my array will change to the new wx value. Did i get the idea correctly? Also my initial conditions are not integer so i need to define an array with float numbers. How can i describe such an array then?

[hknatas]

In here delta_t = time_vals right? And when i tried to run program an error occurs : "Unable to perform assignment because the left and right sides have a different number of elements." What can be the problem?

This is my whole code:

clc;

clear;

clear all

n = 63

Wx(1) = 0.0002 + 0.0001 * n; % 0.0065

Wy(1) = 0.0003 + 0.0001 * n; %0.0066

Wz(1) = 0.0004 + 0.0001 * n; %0.0067

% The initial moments of inertia of the satellite ( m^4 ):

Jx = 2.1e-03;

Jy = 2e-03;

Jz = 1.9e-03;

% The disturbance torque acting on the satellite ( N.m ):

NT = 3.6e-10;

delta_t = 0:0.1:5400;

N_times = length(delta_t);

for i = 1 : N_times

      Wx(i+1) = Wx(i) + delta_t / Jx * (Wz(i) * Wy(i) + NT)*(Jy - Jz);
      Wy(i+1) = Wy(i) + delta_t / Jy * (Wx(i) * Wz(i) + NT)*(Jz - Jx);
      Wz(i+1) = Wz(i) + delta_t / Jz * (Wx(i) * Wy(i) + NT)*(Jx - Jy);

end

Answer 1

[Bruno Luong]

Actually i thought it is going to be easy but when i write the code for wx, i realized that wx is increasing with time but wz and wy in the equation are constant values.

Just copy the value in 3 different variables, then do the calculation.

[hknatas]

I am a beginner. Unfortunately i didn't understand what did you mean by saying copy the calculation. I am just newly learning loops. I should write such a code that: in step one t = 0, program calculates all wx0, wy0 and wz0 values, in step two t = 0.1, program calculates all wx1, wy1 and wz1 values. But in this step program should use the values wx0, wy0 and wz0. And thus, i need to reach 5400s with 54000 iteration. For example in step 3765 program should calculate wx3764 with using wy3763 and wz3763 values. I couldn't figure out which loop can provide this.

[Walter Roberson]

wx(1) = 0.0065
wy(1) = 0.0066
wz(1) = 0.0067
time_vals = 0:0.1:5400;
N_times = length(time_vals);
for i = 1 : N_times
  wz(i+1) = wz(i) + delta_t / Jz * (wx(i) * wy(i) + NT)*(Jx - Jy);
  etc
end

[hknatas]

In here delta_t = time_vals right? And when i tried to run program an error occurs : "Unable to perform assignment because the left and right sides have a different number of elements." What can be the problem?

This is my whole code:

clc;

clear;

clear all

n = 63

Wx(1) = 0.0002 + 0.0001 * n; % 0.0065

Wy(1) = 0.0003 + 0.0001 * n; %0.0066

Wz(1) = 0.0004 + 0.0001 * n; %0.0067

% The initial moments of inertia of the satellite ( m^4 ):

Jx = 2.1e-03;

Jy = 2e-03;

Jz = 1.9e-03;

% The disturbance torque acting on the satellite ( N.m ):

NT = 3.6e-10;

delta_t = 0:0.1:5400;

N_times = length(delta_t);

for i = 1 : N_times

      Wx(i+1) = Wx(i) + delta_t / Jx * (Wz(i) * Wy(i) + NT)*(Jy - Jz);
      Wy(i+1) = Wy(i) + delta_t / Jy * (Wx(i) * Wz(i) + NT)*(Jz - Jx);
      Wz(i+1) = Wz(i) + delta_t / Jz * (Wx(i) * Wy(i) + NT)*(Jx - Jy);

end

[Walter Roberson]

delta_t = 0.1;

time_vals = 0 : delta_t : 5400;
N_times = length(time_vals);

n = 63;

Wx = zeros(1, N_times+1);
Wy = zeros(1, N_times+1);
Wz = zeros(1, N_times+1);
Wx(1) = 0.0002 + 0.0001 * n; % 0.0065
Wy(1) = 0.0003 + 0.0001 * n; %0.0066
Wz(1) = 0.0004 + 0.0001 * n; %0.0067

% The initial moments of inertia of the satellite ( m^4 ):

Jx = 2.1e-03;

Jy = 2e-03;

Jz = 1.9e-03;

% The disturbance torque acting on the satellite ( N.m ):

NT = 3.6e-10;

for i = 1 : N_times

      Wx(i+1) = Wx(i) + delta_t / Jx * (Wz(i) * Wy(i) + NT)*(Jy - Jz);
      Wy(i+1) = Wy(i) + delta_t / Jy * (Wx(i) * Wz(i) + NT)*(Jz - Jx);
      Wz(i+1) = Wz(i) + delta_t / Jz * (Wx(i) * Wy(i) + NT)*(Jx - Jy);

end

plot3(Wx, Wy, Wz);

[hknatas]

I am so grateful, thanks for everthing. One more thing i want to ask: What if i want to Wx as a dependent of time? When i write plot(time_vals,Wx) it says Vectors must be the same length.

[Walter Roberson]

delta_t = 0.1;

time_vals = 0 : delta_t : 5400;
N_times = length(time_vals);

n = 63;

Wx = zeros(1, N_times);
Wy = zeros(1, N_times);
Wz = zeros(1, N_times);
Wx(1) = 0.0002 + 0.0001 * n; % 0.0065
Wy(1) = 0.0003 + 0.0001 * n; %0.0066
Wz(1) = 0.0004 + 0.0001 * n; %0.0067

% The initial moments of inertia of the satellite ( m^4 ):

Jx = 2.1e-03;

Jy = 2e-03;

Jz = 1.9e-03;

% The disturbance torque acting on the satellite ( N.m ):

NT = 3.6e-10;

for i = 1 : N_times - 1

        Wx(i+1) = Wx(i) + delta_t / Jx * (Wz(i) * Wy(i) + NT)*(Jy - Jz);
        Wy(i+1) = Wy(i) + delta_t / Jy * (Wx(i) * Wz(i) + NT)*(Jz - Jx);
        Wz(i+1) = Wz(i) + delta_t / Jz * (Wx(i) * Wy(i) + NT)*(Jx - Jy);

end

plot(time_vals, Wx);

[hknatas]

Thanks for everything, especially for your patience. This really helps a lot.