Solving a 3 variables ODE using a for loop and a step-by-step process

Question

VincentLec 2016-7-27

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/297418-solving-a-3-variables-ode-using-a-for-loop-and-a-step-by-step-process

评论： Stephen23 2016-7-28

Hello,

So here's my problem: I have a 3 dimensional ODE that I need to solve in order to obtain a 3D vector (see attached pdf file for more details). Now I have this for loop relative to the time which allows me to obtains certain components of my ODE. Now, in the end, the c0, c0p, c0pp, psi, psip and psipp (shown in code below) allow me to determine the values of some variables that are used in my ODE. My ODE is theorically a second order ODE by can be considered as a first order since the only unknown variables in the ODE are theta(dotdot) and theta(dot). Also, I defined theta(dotdot) as being thetadotdot(j)=(thetadot(j)-thetadot(j-1))/n where n is the step size. Now, since I wrote that thetadot(1)=[0 ; 0 ; 0], I would've assumed that the command solve would've been enough to isolate thetadot(j) for each time (or each j). However, I always get the message "Undefined variable thetap" which is kinda normal considering that it is the variable that i'm trying to find.

So anyway, do you guys have any idea what I should use to solve my ODE. I already tried doing it all algebraically, but it only works if at least two of the three components of the vector psi are equal to zero (in my case, I want to make it work when all three components of psi are not equal to zero).

Here's part of my code. I've also attached the entire code if it helps.

Thanks a bunch

n=0.01; % Step size
tmax=5; % Maximum time 
for j=2:tmax/n+1  % I used 2 because at j=1, it's time 0 where everything equals 0
  time=n*(j-1); % Time component
  t(:,j)=n*(j-1); % Vector time
  % Translations
  c0(:,1)=[0.05*sin(pi/2*0) ; 0.00*sin(2*pi*0); 0.04*sin(2*pi*0)];             % At time 0
  c0(:,j)=[0.05*sin(pi/2*time) ; 0.00*sin(2*pi*time); 0.04*sin(2*pi*time)];    % MOVEMENT DEFINITION
  c0p(:,j)=(c0(:,j)-c0(:,j-1))/n;                                              % Derivative of position = speed
  c0pp(:,j)=(c0p(:,j)-c0p(:,j-1))/n;                                           % Derivative of speed = acceleration
% Rotations
  psi(:,1)=[0*pi*sin(pi*0) ; 0*pi*sin(pi/2*0) ; 1*pi*sin(pi/4*0)];             % At time 0
  psi(:,j)=[0*pi*sin(pi*time) ; 0*pi*sin(pi/2*time) ; 1*pi*sin(pi/4*time)];    % ANGLES DEFINITION
  psip(:,j)=(psi(:,j)-psi(:,j-1))/n;                                           % Derivative of angular position = angular speed
  psipp(:,j)=(psip(:,j)-psip(:,j-1))/n;                                        % Derivative of angular speed = angular accel.
  ...
  ...
  ...
  ...
  thetap(:,1)=zeros(3,tmax/n+1);
  thp(:,j)=solve((R*(thetap(:,j)-thetap(:,j-1))/n(-tau(:,j)+(I0+I1)*omega0p(:,j)-cross(omega0(:,j),(I0+I1)*omega0(:,j))-(Irx*(omega0p(:,j)+(transpose((thetap(:,j)-thetap(:,j-1))/n)*ex)*(Q*ex)+cross((transpose(thetap(:,j))*ex)*omega0(:,j),Q*ex))+cross((omega0(:,j)+(transpose(thetap(:,j))*ex)*(Q*ex)),Irx*(omega0(:,j)+(transpose(thetap(:,j))*ex)*(Q*ex)))+Iry*(omega0p(:,j)+(transpose((thetap(:,j)-thetap(:,j-1))/n)*ey)*(Q*ey)+cross((transpose(thetap(:,j))*ey)*omega0(:,j),Q*ey))+cross((omega0(:,j)+(transpose(thetap(:,j))*ey)*(Q*ey)),Iry*(omega0(:,j)+(transpose(thetap(:,j))*ey)*(Q*ey)))+Irz*(omega0p(:,j)+(transpose((thetap(:,j)-thetap(:,j-1))/n)*ez)*(Q*ez)+cross((transpose(thetap(:,j))*ez)*omega0(:,j),Q*ez))+cross((omega0(:,j)+(transpose(thetap(:,j))*ez)*(Q*ez)),Irz*(omega0(:,j)+(transpose(thetap(:,j))*ez)*(Q*ez)))))),thetap)
end