Finding the point before landing

1 回表示 (過去 30 日間)
Finnien Mitchell
Finnien Mitchell 2020 年 5 月 24 日
コメント済み: Finnien Mitchell 2020 年 5 月 24 日
I am using a numerical method for finite differences so I can plot the 'x' and 'y' displacement of a projectile neglecting the force of drag. However, as the solutions computed are discrete I can't find the exact point where the projectile lands at (y = 0). So I'm attempting to use the point before/after landing using the code below:
point_before_landing = max(find(fNWy<0));
point_after_landing = max(find(fNWy<0)) + 1;
I am then typing fNWy(point_before_landing) into the command window to find the point. However, when I do this it comes back with the error:
>> fNWy(point_before_landing)
Undefined function or variable 'point_before_landing'.
Would anyone know why this is happening since I believe I defined that variable in my editor section. For reference this is my entire code below as well which is running fine. Any help would be greatly appreciated.
close all
clc
%% Problem Parameters
m = 50*10^-3; % mass in [kg]
radius = 30*10^-3; % radius [m]
a0 = 25; % initial angle [deg]
v0 = 50; % initial velocity [m/s]
y0 = 50; % initial elevation [m]
vy0 = v0*sind(a0); % initial y velocity;
x0 = 0; %initial horizontal
vx0 = v0*cosd(a0); % inital x velocity;
Cdrag = 0.8; % coefficient of drag
Area = 0.0028274334; %cross section area [m^2]
rho = 1.183; % density of air
g = 9.81; % gravity
t = 6.02007; % simulation time [s]
dt = 0.01; % timestep [s]
Niter = floor(t/dt); %iteration
%% Numerical Without Drag (NW)
% To neglect drag Cdrag is set to be equal to 0
Cdrag = 0.00;
NWx(1) = x0;
NWx(2) = x0 + dt*vx0;
NWy(1) = y0;
NWy(2) = y0 + dt*vy0;
NWVx(1) = vx0;
NWVx(2) = vx0;
for i = 3:Niter
NWxNow = NWx(i - 1);
NWxPre = NWx(i - 2);
NWyNow = NWy(i - 1);
NWyPre = NWy(i - 2);
fNWx = @(NWxNex)m*((NWxNex-2*NWxNow+NWxPre)/(dt^2))+((Cdrag*rho*Area)/2) ...
* sqrt(((NWxNow-NWxPre)/dt)^2+ ((NWyNow-NWyPre)/dt)^2) * ...
((NWxNex-NWxPre)/(2*dt));
NWxNex = fzero(fNWx,NWxNow);
NWx(i) = NWxNex;
fNWy = @(NWyNex)m*((NWyNex-2*NWyNow+NWyPre)/(dt^2))+((Cdrag*rho*Area)/2)...
* sqrt(((NWxNow-NWxPre)/dt)^2 +((NWyNow-NWyPre)/dt)^2) * ...
((NWyNex-NWyPre)/(2*dt))+ m*g ;
NWyNex =fzero(fNWy,NWyNow);
NWy(i) = NWyNex;
end
point_before_landing = max(find(fNWy<0));
point_after_landing = max(find(fNWy<0)) + 1;

回答 (1 件)

Thiago Henrique Gomes Lobato
Thiago Henrique Gomes Lobato 2020 年 5 月 24 日
You used the function, not the calculated variables to find the index. Your point definiton was also a bit off. This version of your code works:
close all
clc
%% Problem Parameters
m = 50*10^-3; % mass in [kg]
radius = 30*10^-3; % radius [m]
a0 = 25; % initial angle [deg]
v0 = 50; % initial velocity [m/s]
y0 = 50; % initial elevation [m]
vy0 = v0*sind(a0); % initial y velocity;
x0 = 0; %initial horizontal
vx0 = v0*cosd(a0); % inital x velocity;
Cdrag = 0.8; % coefficient of drag
Area = 0.0028274334; %cross section area [m^2]
rho = 1.183; % density of air
g = 9.81; % gravity
t = 6.02007; % simulation time [s]
dt = 0.01; % timestep [s]
Niter = floor(t/dt); %iteration
%% Numerical Without Drag (NW)
% To neglect drag Cdrag is set to be equal to 0
Cdrag = 0.00;
NWx(1) = x0;
NWx(2) = x0 + dt*vx0;
NWy(1) = y0;
NWy(2) = y0 + dt*vy0;
NWVx(1) = vx0;
NWVx(2) = vx0;
for i = 3:Niter
NWxNow = NWx(i - 1);
NWxPre = NWx(i - 2);
NWyNow = NWy(i - 1);
NWyPre = NWy(i - 2);
fNWx = @(NWxNex)m*((NWxNex-2*NWxNow+NWxPre)/(dt^2))+((Cdrag*rho*Area)/2) ...
* sqrt(((NWxNow-NWxPre)/dt)^2+ ((NWyNow-NWyPre)/dt)^2) * ...
((NWxNex-NWxPre)/(2*dt));
NWxNex = fzero(fNWx,NWxNow);
NWx(i) = NWxNex;
fNWy = @(NWyNex)m*((NWyNex-2*NWyNow+NWyPre)/(dt^2))+((Cdrag*rho*Area)/2)...
* sqrt(((NWxNow-NWxPre)/dt)^2 +((NWyNow-NWyPre)/dt)^2) * ...
((NWyNex-NWyPre)/(2*dt))+ m*g ;
NWyNex =fzero(fNWy,NWyNow);
NWy(i) = NWyNex;
end
[~,point_before_landing] = min(NWy(NWy>0));
point_after_landing = point_before_landing+1;
fNWy(point_before_landing)
  1 件のコメント
Finnien Mitchell
Finnien Mitchell 2020 年 5 月 24 日
thankyou

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

カテゴリ

Help Center および File ExchangeGeometric Geodesy についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by