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 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by