eusquare.m

% Euler solution of second order O.D.E.
% damped spring-mass system with applied square wave forcing function
%
%  ODE is:   m*a + c*v + k*x = f
%
% where f is a square wave of specified amplitude from time t1 to t2

clear;

% define function file smdrk4.m (spring-mass-damper Runge-Kutta 4th order)     
%
%   function a = smdrk4(x,v,f,c,k,m)
%   a = (f - c*v - k*x)/m;

% define constants

m = 1;         % mass
c = 1;         % damping
k = 6;         % spring stiffness

t0 = 0;        % initial time
tmax = 40;     % total time
dt = 0.02;    % time step

t1 = 2;        % start time of square wave
t2 = 14;       % end time of square wave
fsquare = 10;  % square wave amplitude

% initial conditions

x = 0;         % position
v = 0;         % velocity
time = t0;

fprintf(['Solution of Spring-Mass-Damper with Square ',...
         'Wave Forcing Function\n']);
fprintf('Using Euler''s method\n');

% start Euler solution

i = 0;
while time < tmax
  time = time + dt;
  ti = time - dt;     % initial time for this time step
  xi = x;             % initial position for this time step

  if (ti >= t1) & ( ti <= t2)       % evaluate forcing function
    f = fsquare;                    % amplitude
  else
    f = 0;
  end

  xnew = x + v*dt;
  vnew = v + smdrk4(x,v,f,c,k,m)*dt;

  x = xnew;
  v = vnew;

  i = i + 1;
  position(i) = x;
  timev(i) = time;
end;

% the one line Matlab solution; note NEW function: smdode45

[time45,results45] = ode45('smdode45',t0,tmax,[0 0]);

% view results

force = zeros(1,length(timev));  % create force vector for viewing
for i = t1/dt:t2/dt
  force(i) = 1; 
end

% plot the curves

plot(timev,position,time45,results45(:,1),'o',timev,force,'-');
title('Euler and ode45 solution of spring-mass-damper system');
legend('Euler','Matlab''s ode45','forcing function');
axis([t0 tmax -1 3]);