Error in ODE45.

1 回表示 (過去 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 件のコメント
Walter Roberson
Walter Roberson 2021 年 2 月 6 日
%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 ExchangeProgramming についてさらに検索

製品


リリース

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by