Error in ODE45.

Hi guys,

I am trying to simulate a falling object with air resistance and wind involved. I am trying to use ODE45 but I get errors and I dont get why. Please help.

The code:

%set variables:
global g CdS rho uvw m
g = 9.81;
m = 1000;
CdS = 3;
rho = 1.225;
uvw = [1;2;0];
ExitPoint = [0;0;-5000];
GSvExit = [30;40;-7];
Ground = -100;

dt = 0.01;                                          % Time steps (s)
tmax = 100;                                         % Maximum simulation time (s)
tspan = [0:dt:tmax];                                % Simulation time (s)
Initial = [ExitPoint(1),GSvExit(1),ExitPoint(2),GSvExit(2),ExitPoint(3),GSvExit(3)]; % Create an array of initial values for ODE45
options = odeset('AbsTol',1e-6,'MaxStep',5);        % Set some options for ODE45
[t,output] = ode45(@odefcn,tspan,Initial,options);  % Set ODE45 odefcn function

output(output(:,5) > -Ground,:) = [];               % Delete all rows below ground level
t = t(1:length(output));                            % Delete all times that correspond with heigth below ground

function [bal] = odefcn(t,Input)                    % Execute equations of motion
    global g CdS m rho uvw                          % Get global variabeles

    factor = (rho*CdS)./(2*m);                      % Factor to multiply with in equations of motion

    bal = zeros(6,1);                               % Create matrix for results of the equations of motion 
    
    bal(1) = Input(2);                              % horizontal velocity dx/dt (m/s) to location            
    bal(3) = Input(4);                              % horizontal velocity dy/dt (m/s) to location            
    bal(5) = Input(6);                              % vertical velocity dz/dt (m/s) to location
             
    bal(2) = -factor.*(Input(2)-uvw(1)).^2;         % horizontal acceleration d2x/dt2 (m/s) to velocity
    bal(4) = -factor.*(Input(4)-uvw(2)).^2;         % horizontal acceleration d2y/dt2 (m/s) to velocity
    bal(6) = g-factor.*(Input(6)-uvw(3)).^2;        % vertical acceleration d2z/dt2 (m/s) to velocity
end

Answer 1

[Alan Stevens]

You have missed the m in your first global delaration. You should have

global g CdS m rho uvw

not

global g CdS rho uvw

[Walter Roberson]

Right. Or better yet, Don't use global!

http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

[Mike Mierlo van]

Alan and Walter,

I added the global m, but it still does not work.

@Walter. I don’t know how to implement your sugggestion. Is it like: odefcn(t,Input,m,g,Cds,rho,uvw)? I tried it with a working ODE45 like that, but it seems I cannot just add varaiables to odefcn() between the brackets.

So, I don’t have a solution yet.

With all global variables, the error still says:

Index exceeds the number of array elements (0).

Error in Sim>odefcn (line 76) bal(2) = -factor.*(Input(2)-uvw(1)).^2; % horizontal acceleration d2x/dt2 (m/s) to velocity

Regards, Mike

[Mike Mierlo van]

Nevermind... I found another crack in my code.

You did a good job. It is working now.

Thanks a lot!

[Walter Roberson]

%set variables:
g = 9.81;
m = 1000;
CdS = 3;
rho = 1.225;
uvw = [1;2;0];
ExitPoint = [0;0;-5000];
GSvExit = [30;40;-7];
Ground = -100;

dt = 0.01;                                          % Time steps (s)
tmax = 100;                                         % Maximum simulation time (s)
tspan = [0:dt:tmax];                                % Simulation time (s)
Initial = [ExitPoint(1),GSvExit(1),ExitPoint(2),GSvExit(2),ExitPoint(3),GSvExit(3)]; % Create an array of initial values for ODE45
options = odeset('AbsTol',1e-6,'MaxStep',5);        % Set some options for ODE45
[t,output] = ode45(@(t,Input)odefcn(t,Input,g,CdS,rho,uvw,m),tspan,Initial,options);  % Set ODE45 odefcn function

output(output(:,5) > -Ground,:) = [];               % Delete all rows below ground level
t = t(1:length(output));                            % Delete all times that correspond with heigth below ground
plot(t, output); legend({'dx', 'dy', 'dz', 'd2x', 'd2y', 'd2z'})

function [bal] = odefcn(t,Input,g,CdS,rho,uvw,m)                    % Execute equations of motion

    factor = (rho*CdS)./(2*m);                      % Factor to multiply with in equations of motion

    bal = zeros(6,1);                               % Create matrix for results of the equations of motion 
    
    bal(1) = Input(2);                              % horizontal velocity dx/dt (m/s) to location            
    bal(3) = Input(4);                              % horizontal velocity dy/dt (m/s) to location            
    bal(5) = Input(6);                              % vertical velocity dz/dt (m/s) to location
             
    bal(2) = -factor.*(Input(2)-uvw(1)).^2;         % horizontal acceleration d2x/dt2 (m/s) to velocity
    bal(4) = -factor.*(Input(4)-uvw(2)).^2;         % horizontal acceleration d2y/dt2 (m/s) to velocity
    bal(6) = g-factor.*(Input(6)-uvw(3)).^2;        % vertical acceleration d2z/dt2 (m/s) to velocity
end

[Mike Mierlo van]

Thank you Walter!