Error in ODE45.

2 次查看（过去 30 天）
Mike Mierlo van
Mike Mierlo van 2021-2-6
评论： Mike Mierlo van 2021-2-6
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

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

采纳的回答

Alan Stevens
Alan Stevens 2021-2-6
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 个评论
显示 3更早的评论
Walter Roberson
Walter Roberson 2021-2-6
Ran in:
%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 2021-2-6
Thank you Walter!

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

产品


版本

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by