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