Euler’s method
Below is a small script written to demonstrate Euler's method being applied. It simulates a projectile being fired.
The projectile has a known mass, a known drag coefficient, a known initial velocity and a known launch angle. The projectile then expeirences free flight with only drag acting on it. The script terminates when the projectile hits the ground.
Code:
% This is an implementation of Euler's method for solving a differential % problem. It plots the flight path of a projectile experiencing drag. % Inputs: Launch angle, launch velocity % Output: Plot of trajectory % You'll notice that all data is kept, even though it isn't needed. This % has the benefit of making debugging easier, but increases memory usage % and slows the simulation slightly. Best practices suggest that all data % of a simplified simulation should be kept whilst debugging, before % removing storage and increasing complexity of the simulation during % runtime. % CONSTANTS g = 9.81; % Accel. due to gravity, m/s^2 rho = 1.2; % Density of air, kg/m^3 Cd = 1; % Drag coefficient of body dt = 0.01; % Time step, s m = 0.5; % Mass of projectile, kg A = 0.08; % Cross-sectional area, m^2 x = [0]; % X and Y displacements y = [0.001]; % Not 0 so as to begin properly - makes little difference vx = [0]; % X and Y velocities vy = [0]; ax = [0]; % X and Y accelerations ay = [0]; step = 0; % For matrix formation % INPUTS angle = input('Launch angle, degrees ')*pi/180; % Launch angle, rads vel = input('Launch velocity, m/s '); % Initial launch velocity, m/s vx(step+1) = vel*cos(angle); % Velocity in x, m/s vy(step+1) = vel*sin(angle); % Velocity in y, m/s % EULER'S METHOD SOLVER while y > 0 step = step+1; % Step through one time step angle = atan(vy(step)/vx(step)); % calc new angle off velocities vel = sqrt(vx(step)^2+vy(step)^2); % Calc total vel. for drag. m/s drag = vel*Cd*A; % Find total drag force, N ax(step+1) = -1*drag*cos(angle)/m; % Accel. in x, m/s^2 ay(step+1) = -1*drag*sin(angle)/m - g; % Accel. in y, m/s^2 vx(step+1) = vx(step) + ax(step+1)*dt; % Vel. update, m/s vy(step+1) = vy(step) + ay(step+1)*dt; % Vel. update, m/s x(step+1) = x(step) + vx(step+1)*dt + 0.5*ax(step+1)*dt^2; % Pos. update, m y(step+1) = y(step) + vy(step+1)*dt + 0.5*vy(step+1)*dt^2; % Pos. update, m end plot(x,y); % Plot trajectory