rksquare.m

% Runge-Kutta 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.1;      % 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 manual 4th order Runge-Kutta and Matlab''s ode45\n');

% start Runge-Kutta 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

  k1x = v;
  k1v = smdrk4(x, v,f,c,k,m);

  k2x = v + k1v * dt/2;
  k2v = smdrk4(x + k1x * dt/2, v + k1v * dt/2,f,c,k,m);

  k3x = v + k2v * dt/2;
  k3v = smdrk4(x + k2x * dt/2, v + k2v * dt/2,f,c,k,m);

  k4x = v + k3v * dt;
  k4v = smdrk4(x + k3x * dt, v + k3v * dt,f,c,k,m);

  x = x + (k1x + 2*k2x + 2*k3x + k4x) * dt/6;
  v = v + (k1v + 2*k2v + 2*k3v + k4v) * dt/6;

  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('Runge-Kutta solution of spring-mass-damper system');
legend('4th order Runge-Kutta','Matlab''s ode45','forcing function');
axis([t0 tmax -1 3]);