Error in ODE45.

10 views (last 30 days)
Mike Mierlo van
Mike Mierlo van on 6 Feb 2021
Commented: Mike Mierlo van on 6 Feb 2021
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

Accepted Answer

Alan Stevens
Alan Stevens on 6 Feb 2021
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
  5 Comments
Walter Roberson
Walter Roberson on 6 Feb 2021
%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
Mike Mierlo van on 6 Feb 2021
Thank you Walter!

Sign in to comment.

More Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2020b

Community Treasure Hunt

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

Start Hunting!