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

Specific points to note are as follows:

1. Adjusting the time step will very quickly make it inaccurate. Remember that it effectively simulates a lot of little straight displacements, and this is demonstrated if you make the time step around 0.1 seconds, using a 30 degree angle and a 10m/s initial velocity. Because of this, it's very brute-force and requires a lot of calculation cycles.

2. Note that you can easily recognise the formulas to calculate velocities and displacements. v = u +a*t, and s = s[sub]0[/sub] + u*t + 0.5*a*t^2. These are just the linear motion formulas we're all familiar with.

It all works off finding integrals, as explained in the parent document. In this case, we exploit the relationship that velocity and acceleration of an object are the first and second derivatives of its displacement, respectively, and work backwards, as we can solve for acceleration directly. However, remember that this is not the only capability of the method indeed, it can solve many sorts of integral problems.

Possible Related Links

• Firearms