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

7 views (last 30 days)
VincentLec on 27 Jul 2016
Commented: Stephen Cobeldick on 28 Jul 2016
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

#### 1 Comment

Stephen Cobeldick on 28 Jul 2016

Torsten on 28 Jul 2016
I wonder why you don't use one of MATLAB's ODE-integrators, e.g. ODE15S.
Best wishes
Torsten.