Input result derivative into ode45

2020 4 25
1 个回答
回答已采纳
更新时间：2020 4 26
6 次查看（30 天）

0 个投票

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)

6 个评论
显示 4更早的评论 隐藏 4更早的评论

darova
darova 2020-4-25
Is it possible to see formulas?
Marc
Marc 2020-4-25
编辑：Marc 2020-4-25
%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
darova 2020-4-25
Can you show original formulas?
like
Marc
Marc 2020-4-25
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
darova 2020-4-25
You have , , , , , (6 uknowns)
But i only see 4 equations. Do you have more?
Marc
Marc 2020-4-25
According to the problem modelling, y and x are constant, and bother their derivatives and second derivatives are zero.

请先登录，再进行评论。

请先登录，再回答此问题。

 采纳的回答

darova
darova 2020-4-25

0 个投票

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

3 个评论
显示 1更早的评论 隐藏 1更早的评论

Marc
Marc 2020-4-26
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
Marc 2020-4-26
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
darova 2020-4-26
Yes you understood me correctly

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by