Using finite differences to find the displacement (x,y) of a projectile including the effect of drag
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponses (1)
KSSV
le 22 Mai 2020
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
7 commentaires
KSSV
le 23 Mai 2020
The code is working fine for me......Breakingpoint can be easily known...it will appear as a red dot on the left of the editor.
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!