function ex11
%  This is a mechanics problem studied in section 4 of Chapter XII of
%  R. Mattheij and J. Molenaar, Ordinary Differential Equations in Theory
%  and Practice, Wiley, New York, 1996.  It describes the motion of a bead 
%  on a rotating needle subject to the forces of gravity, friction, and 
%  centrifugal force.  The equations of motion are
%
%    u1' =   u2
%    u2' = - fr*u2 + lambda*psi(t)
%    u3' =   u4
%    u4' = - fr*u4 - lambda + g/omega^2
%
%  There is a kinematic constraint that corresponds to the bead having
%  to move along the needle:  0 = u3 - u1*psi(t)
%
%  Here the coordinates of the bead are (u1,u3).  The quantity lambda
%  is a Lagrange multiplier.  The function psi(t) = tan(t + phi0), 
%  where phi0 is the initial angle between the needle and the horizontal
%  that is taken here to be pi/4.  The needle is rotating with a constant 
%  angular velocity omega and g is the acceleration due to gravity.  Here
%  this ratio is taken to be 1.  fr is a non-dimensional friction parameter
%  that is taken here to be 10.  The bead is started from (u1,u3) = (1,1) 
%  with initial derivatives (u2,u4) = (-6,-6).
%
%  It is inconvenient that psi becomes infinite during the integration,
%  so we introduce the variable u5 = lambda/cos(t+phi0) in the equations
%  of motion and rewrite the kinematic constraint to obtain the DAE
%
%    u1' =   u2
%    u2' = - 10*u2 + u5*sin(t+pi/4)
%    u3' =   u4
%    u4' = - 10*u4 - u5*cos(t+pi/4) + 1
%     0  =   cos(t+pi/4)*u3 - sin(t+pi/4)*u1 = g(t)
%
%  Replacing the kinematic constraint 0 = g(t) by the implicit constraint
%  0 = g''(t) results in a DAE of index 1.  Unfortunately, as is common with
%  mechanics problems, when this is done the numerical solution drifts so 
%  that the constraint 0 = g(t) is not satisfied well and the computation 
%  breaks down. Here Baumgarte stabilization is done as explained on p. 470 
%  of E. Hairer and G. Wanner, Solving Ordinary Differential Equations II
%  Stiff and Differential-Algebraic Problems, 2nd ed., Springer, Berlin,
%  1996. In this approach the equation is replaced not by 0 = g''(t) but by
%
%    0 = g''(t) + 2*alpha*g'(t) + alpha^2 *g(t)
%
%  for a suitable value of alpha > 0.  Here alpha is taken to be 10 and a 
%  little calculation shows that
%
%    g'(t)  =   cos(t+pi/4)*(u4 - u1) + sin(t+pi/4)*(- u2 - u3)
%    g''(t) =   cos(t+pi/4)*(u4' - u1' - u2 - u3) 
%             + sin(t+pi/4)*(- u2' - u3' - u4 + u1)

% Copyright 2002, The MathWorks, Inc.

u0 = zeros(5,1);
u0(1) =  1;
u0(2) = -6;
u0(3) =  1;
u0(4) = -6;
u0(5) = -10;       %  A guess for the algebraic variable.

M = diag([1 1 1 1 0]);
options = odeset('Mass',M,'RelTol',1e-4);

[t,u] = ode15s(@ex11ode,[0 15],u0,options);

plot(u(1,1),u(1,3),'o',u(:,1),u(:,3))
title('Path of a bead on a rotating needle.')

%-------------------------------------------------------------------------%

function up = ex11ode(t,u)

  s = sin(t+pi/4);
  c = cos(t+pi/4);
  up = zeros(5,1);
  up(1) =   u(2);
  up(2) = - 10*u(2) + u(5)*s;
  up(3) =   u(4);
  up(4) = - 10*u(4) - u(5)*c + 1;

  g     =  c*u(3) - s*u(1);
  gp    =  c*(u(4) - u(1)) + s*( - u(2) - u(3));
  gpp   =  c*(up(4) - up(1) - u(2) - u(3)) ...
         + s*(- up(2) - up(3) - u(4) + u(1));
  up(5) = gpp + 20*gp + 100*g;
    