Solving a differential matrix ode with ode45

I am trying to estimate the angular velocity of a satellite given some input torques, its previous angular velocity, and other characteristics of the satellite. I'm having problems passing in the previous angular velocity vector (the code that isn't working is attached).

here is the function describing the dynamics:

function w_body_dot = dynamics(w_body, flag, P)
% model of the satellite dynamics
% P is [J_satellite; Ai; J_w; Wheel_velocities; input_torques] with two
% trailing columns of zeros on wheel velocities and input torques
J_satellite = P(1:3,1:3);
Ai = P(4:6,1:3);
J_w = P(7:9,1:3);
Wheel_velocities = P(10:12,1);
input_torques = P(13:15,1);
J_hat = J_satellite - Ai*J_w*Ai';
J_hat_mult_w_body_dot = cross(-w_body, (J_satellite*w_body + Ai*J_w*Wheel_velocities)) + Ai*input_torques;
w_body_dot = linsolve(J_hat, J_hat_mult_w_body_dot);
end

P is used to pass in values that I'm not solving for, but I want to be able to change and determine the response. let's have:

P =  [0.0604         0         0;
           0    0.0420         0;
           0         0    0.0880;

      1.0000         0         0;
           0    1.0000         0;
           0         0    1.0000;

      0.0000         0         0;
           0    0.0000         0;
           0         0    0.0000;

           0         0         0;
           0         0         0;
           0         0         0;

           0         0         0;
           0         0         0;
           0         0         0]

Then I choose a starting w_body, some time step, and try to run ode45:

w_body = [0;0;0];
sample_time = 0.1
[t, w_body] = ode45(@dynamics, [0 sample_time], w_body, [], P);

I get the following error:

Error using  + 
Matrix dimensions must agree.

I found this was because as soon as I pass w_body into ode45, it gets cut to be just 0. This causes problems for the matrices in the dynamics function. How can I pass in a vector value for w_body?

I tried passing in w_body as a horizontal vector instead, and then transposing it in the dynamics function, but this didn't work (ode45 still turned it into just 0).

[Pinky]

I realized I had left a line out of the dynamics function ( to calculate J_hat), and have updated the question to include this line.

Answer 1

[Francisco J. Triveno Vargas]

HI my friend, i did two or three corrections now is running. You need to check if has sense or not.

Regards

Francisco

function Test()

%UNTITLED4 Summary of this function goes here

% Detailed explanation goes here

P = [0.0604 0 0;

0 0.0420 0;

0 0 0.0880;

1.0000 0 0;

0 1.0000 0;

0 0 1.0000;

0.0000 0 0;

0 0.0000 0;

0 0 0.0000;

0 0 0;

0 0 0;

0 0 0;

0 0 0;

0 0 0;

0 0 0];

w_body_in = [0.1;0;0.1]; %initial condition

sample_time = 0.01;

[t, w_body] = ode45(@dynamics, [0:sample_time:5], w_body_in, [], P);

figure

plot(t,w_body)

grid on

end

function w_body_dot = dynamics(flag,w_body, P)

% model of the satellite dynamics

% P is [J_satellite; Ai; J_w; Wheel_velocities; input_torques] with two

% trailing columns of zeros on wheel velocities and input torques

J_satellite = P(1:3,1:3);

Ai = P(4:6,1:3);

J_w = P(7:9,1:3);

Wheel_velocities = P(10:12,1);

input_torques = P(13:15,1);

w_body;

J_hat = J_satellite - Ai*J_w*Ai';

J_hat_mult_w_body_dot = cross(-w_body, (J_satellite*w_body + Ai*J_w*Wheel_velocities)) + Ai*input_torques;

w_body_dot = linsolve(J_hat, J_hat_mult_w_body_dot);

end