# Error in ODE45.

1 ビュー (過去 30 日間)
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 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 件のコメント表示非表示 4 件の古いコメント
Mike Mierlo van 2021 年 2 月 6 日
Thank you Walter!

サインインしてコメントする。

R2020b

### Community Treasure Hunt

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

Start Hunting!

Translated by