Input result derivative into ode45

Good morning,

I am trying to implement a code that given angular velocity components, its derivative and 3 initial angles integrates the value of the angular velocity using ode45. The problem is that there is variable dissipation involved, and to calculate the dissipation in each step of time you require both the angular velocity and its derivative in that moment of time. I have tried gradient() but since the step size is variable it doesnt work. The other option is to calculate the angles at each point of time and use them to generate the mass matrix and find the derivative value, but since it is a 2-mass system i cannot find the equations for that (only for single mass systems)

[darova]

Is it possible to see formulas?

[Marc]

%Calculate w_dot0 given the angles in quaternions and initial w
for i = 1:length(t)
Q = [q(i,1)^2-q(i,2)^2-q(i,3)^2+q(i,4)^2,     2*(q(i,1)*q(i,2)+q(i,3)*q(i,4)),     2*(q(i,1)*q(i,3)-q(i,2)*q(i,4));
   2*(q(i,1)*q(i,2)-q(i,3)*q(i,4)),        -q(i,1)^2+q(i,2)^2-q(i,3)^2+q(i,4)^2, 2*(q(i,2)*q(i,3)+q(i,1)*q(i,4));
   2*(q(i,1)*q(i,3)+q(i,2)*q(i,4)),         2*(q(i,2)*q(i,3)-q(i,1)*q(i,4)),    -q(i,1)^2-q(i,2)^2+q(i,3)^2+q(i,4)^2
  ];

M = Q*[-Mass*g*d*Q(3,2);
        Mass*g*d*Q(3,1);
        0];

wx_dot0(i) = M(1,1)/Ix-(Iz-Iy)*wy(i)*wz(i)/Ix;
wy_dot0(i) = M(2,1)/Iy-(Ix-Iz)*wz(i)*wx(i)/Iy;

F0(i,1) = a*wy_dot0(i)-b*wx_dot0(i)-a*wz(i)*wx(i)-b*wy(i)*wz(i);
end

%Fourier fit of the data
F = fit(t,F0,'fourier5');
cf = coeffvalues(F);

%Fourier An and Bn
An = zeros(1,(length(cf)-2)/2);
Bn = zeros(1,(length(cf)-2)/2);
for i = 1:(length(cf)-2)/2
    An(i) = cf(2*i);
    Bn(i) = cf(2*i+1);
end

% Terms for finding z and derivatives
for i = 1:length(An)
    Dn(i) = (An(i)^2+Bn(i)^2)^(1/2);
    Sn(i) = i*cf(end) ;
    En(i) = Dn(i)/(((c2-Sn(i)^2)^2+(c1*Sn(i))^2)^(1/2));
    phin(i) = atan(An(i)/Bn(i))-atan((c1*Sn(i)/(c2-Sn(i)^2)));
end

for i = 1:length(An)
    zn(i) = En(i)*sin(phin(i));
    z_dotn(i) = zn(i)*Sn(i);
end
%z and derivatives for each point
z = sum(zn);
z_dot = sum(z_dotn);
z_dotdot = F0(1)-c1*z_dot-c2*z;  

%Movement equations
A = [Ix+mu*(y^2+z^2),   -mu*x*y,           -mu*x*z;
     -mu*y*x,           Iy+mu*(z^2+x^2),   -mu*z*y;
     -mu*z*x,           -mu*z*y,           Iz+mu*(x^2+y^2)
    ];

B = [-(Iz-Iy+mu*(y^2-z^2))*wy*wz-mu*((2*y*y_dot+2*z*z_dot)*wx+y*z*(wz^2-wy^2)-2*x_dot*y*wy-2*x_dot*z*wz-x*z*wz*wy+x*y*wx*wz+y*z_dotdot-z*y_dotdot);
     -(Ix-Iz+mu*(z^2-x^2))*wz*wx-mu*((2*z*z_dot+2*x*x_dot)*wy+z*x*(wx^2-wz^2)-2*y_dot*z*wz-2*y_dot*x*wx-y*x*wy*wz+y*z*wy*wx+z*x_dotdot-x*z_dotdot);
     -(Iy-Ix+mu*(x^2-y^2))*wx*wy-mu*((2*x*x_dot+2*y*y_dot)*wz+x*y*(wy^2-wx^2)-2*z_dot*x*wx-2*z_dot*y*wy-z*y*wz*wx+z*x*wz*wy+x*y_dotdot-y*x_dotdot);
    ];
sol = A\B;

wx_dot = sol(1)
wy_dot = sol(2)
wz_dot = sol(3)

It is part of a larger code, the 3 w()_dot and z, z_dot, z_dotdot are the solutions i want ode45 to give me, but as it can be seen you need F (and consequently F0) to find the values of z, and f0 requires w_dot. Everything but q, w, w_dot and z, z_dot, z_dotdot is constant

[darova]

Can you show original formulas?

like

[Marc]

Cannot show formulas for the initial 'q' part since those are valid just for the initial condition, after that the formulas change and i haven't been able to find the new ones. However, only the ones starting from the F0(i,1) should be necessary

[darova]

You have , , , , , (6 uknowns)

But i only see 4 equations. Do you have more?

[Marc]

According to the problem modelling, y and x are constant, and bother their derivatives and second derivatives are zero.

Answer 1

[darova]

Here is an idea

Where M X F as following (sorry couldn't finish matrix, it's too long)

Then just put all those matrices inside ode45 function

function du = f(t,u)
wx = u(1);
wy = u(2);
wz = u(3);
z = u(4);
dz = u(5);
% and other parametrs and constants

M = [ ... long matrix ]
F = [ ... long matrix ]
X = M\F;                            % solve matrix
du(1:4,1) = [X(1:3); dz; X(4)];     % pass dwx dwy dwz dz ddz
end

[Marc]

So, you mean letting the first 3 rows of the M and F matrices be the coefficients of the four variables and the independent term extracted from the formulas of the first two pictures and then building the last row with the equations for F below in the same way?

Leaving M like this

M = [Ix+mu*(y^2+z^2),   -mu*x*y,           -mu*x*z          mu*y;
     -mu*y*x,           Iy+mu*(z^2+x^2),   -mu*z*y          -mu*x;
     -mu*z*x,           -mu*z*y,           Iz+mu*(x^2+y^2)  0;
     b,                  -a,               0,               1;   
    ]

Makes sense for me, but i am confused on why the book presents me the rest of equations in order to solve the problem.

[Marc]

Ok, i'm not getting the exact numerical results i expected but the shapes of the envelopes of the curves are right, thank you so much!

[darova]

Yes you understood me correctly