This is an illustration of how the numerical integrator Shake and Rattle work on a simple planar pendulum example. The way it works is by kicking the mass upwards, releasing the constraint, let the mass fly and land back on the constrained set. The drawing is geometrically correct. For more details on those methods, see http://arxiv.org/abs/1207…

gistfile1.tex

\documentclass{standalone}


\usepackage{tikz}
\usetikzlibrary{calc,decorations.markings,intersections,fpu}

\tikzset{
  >=stealth,
  fiber/.style={draw,thick,blue!70!black},
  flight/.style={draw,black,thick, densely dotted, decoration={markings, mark=at position 0.6 with {\arrow{>}}}, postaction={decorate}},
  rattle/.style={color=green!50!black},
}



\begin{document}

  \begin{tikzpicture}[
  >=stealth,
  scale=2.5,
  hinge/.style={fill=white,draw=black},
  fiber_projection/.style={opacity=.3, ->},
  ]


\coordinate (origin) at (0,0);


\newcommand*\pathfreefall[1]{
\path[name path=freefall, #1] (2,-3.5) parabola[parabola height=2cm] +(-3,0);
}
\newcommand*\pathconstraint[1]{
\path[name path=constraint, #1] (0,0) circle[radius=2cm];
}

% Compute the free fall trajectory by clipping a circle and a parabola
\begin{scope}
\pathconstraint{clip}
\pathfreefall{draw, flight}
\end{scope}


% Just "path" the circle and parabola, to compute intersections
\pathconstraint{}
\pathfreefall{}


% Find intersections of the circle and the parabola
\path [name intersections={of=freefall and constraint, by={q0,q1}}];

% The clipping circles are too big: reset the bounding box
\pgfresetboundingbox

% Draw just a bit of the circle which is clipped at a circle centered at the midpoint between q0 and q1
\begin{scope}
	\path[clip] ($(q0)!.5!(q1)$) circle[radius=1cm];
	\pathconstraint{draw}
\end{scope}


% Detla x to compute the tangents
\newcommand*\deltax{.005}
% Vector length
\newcommand*\veclength{.5}

% Compute multiplication factor to obtain a tangent vector of right size
\pgfkeys{/pgf/fpu, /pgf/fpu/output format=fixed}
\pgfmathsetmacro\multfactor{\veclength/\deltax}
\pgfkeys{/pgf/fpu=false}

% Common command for q0 and q1
\newcommand*\drawz[1]{
% draw imaginary line slight off by deltax
\path[name path=vert#1] ($(q#1) + (-\deltax,-1)$) -- ($(q#1) + (-\deltax,1)$);
% compute intersection wiht parabola: this is a base for the derivative
\fill [name intersections={of=freefall and vert#1, by={q#1p}}] ;

% prolong that vector \multfactor times, this gives us the launching/crashing velocity
\coordinate (z#1p) at ($(q#1)!\multfactor!(q#1p)$);
% draw it
\draw[->] (q#1) -- (z#1p);

% orthogonal projection of velocity: this is the reaction force
\coordinate (q#1v) at ($(origin)!(z#1p)!(q#1)$);
% velocity tangent to the circle:
\draw[rattle,thick,->] (q#1) -- +($(z#1p) - (q#1v)$) node[shape=coordinate] (z#1) {};
}

\drawz{0};


% label z_0^+
\node[left] at (z0p) {$z_0^+$};

\drawz{1};

% draw two past and future versions of the rod
\draw (origin) -- (q0);
\draw[gray] (origin) -- (q1);

% label z_1^-
\node[right] at (z1p) {$z_1^-$};
% label z_1
\node[below right] at (z1) {$z_1$};


% create point z which is on \Man but not on \Mp
\coordinate (z) at ($(z0p)!1.3!(z0)$);
% draw z
\draw[->] (q0) -- (z);
% label z
\node[below] at (z) {$z$};

% draw projection from launching/crashing velocity to tangent velocity
\draw[fiber, fiber_projection] (z) -- (z0p);
\draw[fiber, fiber_projection] (z1p) -- (z1);

% draw initial "kick" reaction force sending z to z_0^+
\path[fiber,thick,draw,->] (q0) -- +($(z0p) - (z)$);

% draw final Rattle "kick"
\path[fiber,thick,draw,->] (q1) -- +($(z1) - (z1p)$);


% draw nice hinge
\draw[hinge] (origin) circle(1pt);
\draw[fill=black] (origin) circle(.5pt);


\end{tikzpicture}
\end{document}