Using finite differences to find the displacement (x,y) of a projectile including the effect of drag

Question

Finnien Mitchell 2020-5-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/530148-using-finite-differences-to-find-the-displacement-x-y-of-a-projectile-including-the-effect-of-drag

编辑： darova 2020-5-23

When i enter the code below I get the following error:

Index exceeds the number of array elements (2).
Error in EGB211_CLA (line 37)
    XxNow = XxN(i-1) ;

Would anyone know how to fix this issue? My code is below as well

close all
clc
 
%% Problem Parameters 
 
m = 0.05 ;  % mass (kg) = 50x10^-3
angle = 25 ; % initial angle (degrees)
th0 = angle*(pi/180) ; % initial angle (radians)
Xx0 = 0 ; % initial horizontal displacement (m)
Xy0 = 50 ; % initial vertical displacement (m)
V0 = 50 ; % initial velocity magnitude (m/s)
Vx0 = 41.315 ; % initial horizontal velocity (m/s)
Vy0 = 21.1309 ; % initial vertical velocity (m/s)
r = 0.03 ; % radius (m)
A = pi*(r)^2 ; % cross sectional area of projectile (m^2)
Cdrag = 0.80 ; % drag coefficient (unitless)
rho = 1.183 ; % air density (kg/m^3)
g = 9.81 ; % gravity (m/s^2)
 
 
%% Numerical Solution (With Drag) 
 
t = 10 ;
dt = 0.01 ; % Time step (seconds)
Niter = floor(10/0.01) ; % total number of intervals (t/dt)
 
XxN(1) = Xx0 ; 
XxN(2) = Xx0+(dt*Vx0) ; 
    
XyN(1) = Xy0 ; 
XyN(2) = Xy0+(dt*Vy0) ; 
 
for i = 3:Niter
    
    XxNow = XxN(i-1) ; 
    XxPrev = XxN(i-2) ;
    
    XyNow = XyN(i-1) ; 
    XyPrev = XyN(i-2) ;
    
    fx = @(XxNext)m*((XxNext-(2*XxNow)+XxPrev)/(dt^2)+(0.5*Cdrag*rho*A)...
        *(sqrt(((XxNow-XxPrev/dt)^2)+((XyNow-XyPrev/dt)^2))))*(XxNext-XxPrev/2*dt) ;
    
    fy = @(XyNext)m*((XyNext-(2*XyNow)+XyPrev)/(dt^2)+(0.5*Cdrag*rho*A)...
        *(sqrt(((XxNow-XxPrev/dt)^2)+((XyNow-XxPrev/dt)^2))))*(XyNext-XyPrev/2*dt)+(m*g) ;
    
    XxNext = fzero(fx,XxNow) ; 
    xxN(i) = XxNext ; 
    
    XyNext = fzero(fy,XyNow) ; 
    xyN(i) = XyNext ; 
    
end 

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

KSSV 2020-5-22

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/530148-using-finite-differences-to-find-the-displacement-x-y-of-a-projectile-including-the-effect-of-drag#answer_435993

在 MATLAB Online 中打开

There was a type error in updating XxN.....it was typed as xxN instead of XxN.....find the below corrected code:

close all
clc
%% Problem Parameters 
m = 0.05 ;  % mass (kg) = 50x10^-3
angle = 25 ; % initial angle (degrees)
th0 = angle*(pi/180) ; % initial angle (radians)
Xx0 = 0 ; % initial horizontal displacement (m)
Xy0 = 50 ; % initial vertical displacement (m)
V0 = 50 ; % initial velocity magnitude (m/s)
Vx0 = 41.315 ; % initial horizontal velocity (m/s)
Vy0 = 21.1309 ; % initial vertical velocity (m/s)
r = 0.03 ; % radius (m)
A = pi*(r)^2 ; % cross sectional area of projectile (m^2)
Cdrag = 0.80 ; % drag coefficient (unitless)
rho = 1.183 ; % air density (kg/m^3)
g = 9.81 ; % gravity (m/s^2)
%% Numerical Solution (With Drag) 
t = 10 ;
dt = 0.01 ; % Time step (seconds)
Niter = floor(10/0.01) ; % total number of intervals (t/dt)
XxN = zeros(Niter,1) ; 
XyN = zeros(Niter,1) ; 
XxN(1) = Xx0 ; 
XxN(2) = Xx0+(dt*Vx0) ; 
XyN(1) = Xy0 ; 
XyN(2) = Xy0+(dt*Vy0) ; 
for i = 3:Niter
  
  XxNow = XxN(i-1) ; 
  XxPrev = XxN(i-2) ;
  
  XyNow = XyN(i-1) ; 
  XyPrev = XyN(i-2) ;
  
  fx = @(XxNext)m*((XxNext-(2*XxNow)+XxPrev)/(dt^2)+(0.5*Cdrag*rho*A)...
  *(sqrt(((XxNow-XxPrev/dt)^2)+((XyNow-XyPrev/dt)^2))))*(XxNext-XxPrev/2*dt) ;
  
  fy = @(XyNext)m*((XyNext-(2*XyNow)+XyPrev)/(dt^2)+(0.5*Cdrag*rho*A)...
  *(sqrt(((XxNow-XxPrev/dt)^2)+((XyNow-XxPrev/dt)^2))))*(XyNext-XyPrev/2*dt)+(m*g) ;
  
  XxNext = fzero(fx,XxNow) ; 
  XxN(i) = XxNext ; 
  
  XyNext = fzero(fy,XyNow) ; 
  XyN(i) = XyNext ; 
  
end 