I've been introduced to Mathematica very recently. Basically, I haven't actually "solved" anything in my Mathematica lifetime, but I've done some simulations. With my fractional knowledge, I tried simulating a toggle mechanism today. Here's what I did...

ω0 = .02; ω = .02; r = 2; r1 = 6; r0 = 1.5;
rA[θ_, r_] := r*{Cos[θ], Sin[θ]};
rAy[θ_, r_] := r*Sin[θ];
rAx[θ_, r_] := r*Cos[θ];
rC[θ0_, r0_] := r0*{Cos[θ0] + 1/r0, Sin[θ0] + 3/r0};
rCy[θ0_, r0_] := r0*Sin[θ0] + 3;
rCx[θ0_, r0_] := r0*Cos[θ0] + 1;
xB[θ_, r_] := r*Cos[θ] + Sqrt[(r1^2 - (r*Sin[θ])^2)];
rB[θ_, r_] := {xB[θ, r], 0};
rBx[θ_, r_] := xB[θ, r];
rBy[θ_, r_] := 0;
Manipulate[θ = ω*t; θ0 = ω0*t; 
 rAO = rA[θ, r]; rBO = rB[θ, r]; rDC = rC[θ0, r0];
 Graphics[{{Thick, Darker[Green], Dashed, Circle[{0, 0}, r]}, {Thick, 
    Dashed, Blue, Circle[{1, 3}, r0]},
   {Red, EdgeForm[{Thick, Black}], White, 
    Rectangle[{rBO[[1]] - 0.5, rBO[[2]] - 0.5}, {rBO[[1]] + 1, 
      rBO[[2]] + 0.5}]},
   {Thick, Darker[Red], Line[{{0, 0}, rAO, rBO}]}, {Thick, 
    Darker[Red], Line[{{1, 3}, rDC, rAO}]},
   {Black, Disk[{0, 0}, 0.1], Disk[{1, 3}, 0.1], 
    Disk[{rAO[[1]], rAO[[2]]}, 0.1], Disk[{rDC[[1]], rDC[[2]]}, 0.1], 
    Disk[{rBO[[1]], rBO[[2]]}, 0.1]}},
  PlotRange -> {{-3, 9}, {-3, 5}}, 
  ImageSize -> {600, 250}], {{t, 20, "Motion"}, 20, N[205*Pi], 

I think you can figure out that there's a big issue, due to my enormous expectation from the "visualizing" without even thinking of "solving"...

Both the cranks have the same angular velocity. Because, that's the way I've defined. I've got no other choice. But, the mechanism should begin with the top crank. That crank makes the bottom one to rotate, by transferring power through the connecting rod.

Given a connecting rod of constant length (especially, less than the sum of radii of the cranks), the bottom crank cannot complete one full rotation. It can go only halfway, and then turn back. The slider motion remains the same, but the crank is causing a big trouble to the mechanism due to the ridiculous "elastic behavior" of the connecting rod.

I don't want to analyze this mechanism. So, I'm cool with any answer that essentially simulates the actual mechanism (even if that's a modification of mine) as I just need the simulation only for a video (at least, for now). At the same time, I'll be very happy with "solving" the equations with constraints, because I guess those solutions may help me in the future for other mechanisms.

  • $\begingroup$ I think you first have to choose the length of the rod so that it fits in the distance between the two centers, and subsequently use trigonometry to calculate the position of the point on the lower circle. $\endgroup$ Jan 16, 2014 at 18:22

3 Answers 3


One way is to set up a DAE: See tutorial/DSolveExamplesOfDAEs and example/ModelConstrainedSystemsAsDAEs.

The constraint that the driver (bottom rotating link) has a fixed length is taken care of by initial conditions and the DE. There are two possible starting positions for the driven link. One might have to inspect the result of Solve to determine which is desired.

constraints = ComplexExpand@{
    (* SquaredEuclideanDistance[{a1[t], a2[t]}, {0, 0}] == r1^2, not needed *)
    SquaredEuclideanDistance[{b1[t], b2[t]}, {a1[t], a2[t]}] == rod^2,
    SquaredEuclideanDistance[{c1, c2}, {b1[t], b2[t]}] == r2^2};

eqn = {
      (* DE to rotate small link *)
    D[{a1[t], a2[t]}, t] == omega/r1 Cross[{a1[t], a2[t]}],
      (* initial conditions *)
    {a1[0], a2[0]} == {0, r1},
    {b1[0], b2[0]} == {b1[t], b2[t]} /.
      First @ Solve[constraints /. Thread[{a1[t], a2[t]} -> {0, r1}], {b1[t], b2[t]}]};


Block[{omega = 1, r1 = 1, rod = 3, r2 = 2, c1 = 1, c2 = 3},
  sol = First @ NDSolve[{eqn, constraints},
     {a1[t], a2[t], b1[t], b2[t]},
     {t, 0, 2 Pi/omega}]];

With[{c1 = 1, c2 = 3},
   {Line[{{0, 0}, {a1[t], a2[t]}, {b1[t], b2[t]}, {c1, c2}} /. sol /. t -> t0]},
   PlotRange -> {{-2, 6}, {-2, 6}}
  {t0, 0, Dynamic@sol[[1, 2, 0, 1, 1, 2]]}]

enter image description here

[The piston can be included via trigonometry or by adding another constraint.]

A random linkage:


Module[{sol, x, y, t},
 DynamicModule[{wheels = 6, omega = 1, joints, centers, constraints, 
   variables, initial, eqn, t0 = 0., direction = Forward},

  joints = 
   Sort@RandomReal[1, {wheels, 2}] + Table[{i, 0}, {i, wheels}];
  centers = 
   Sort@RandomReal[1, {wheels, 2}] + Table[{i, (-1)^i}, {i, wheels}];
  variables = 
   Transpose@{Array[x[#][t] &, wheels], Array[y[#][t] &, wheels]};
  constraints = ComplexExpand@{
      SquaredEuclideanDistance[#1, #2] == 
        SquaredEuclideanDistance[#3, #4] &,
      {Most[variables], Rest[variables], Most[joints], Rest[joints]}
      SquaredEuclideanDistance[#1, #3] == 
        SquaredEuclideanDistance[#2, #3] &,
      Rest /@ {variables, joints, centers}
  initial = Thread /@ Thread[variables == joints /. t -> 0] // Flatten;
  eqn = {D[First@variables, t] == 
        First@centers] Cross[First@variables - First@centers]};

  Check[sol = 
    First@NDSolve[Flatten@{eqn, constraints, initial}, 
      Flatten[variables], {t, -2 Pi/omega, 2 Pi/omega}], 
   direction = ForwardBackward, {NDSolve::ndsz, NDSolve::ndcf}];

     Manipulator[Dynamic[t0], First[x[1][t] /. sol /. t -> "Domain"], 
      AnimationDirection -> direction], "t", Left],
     Dynamic@With[{pts = variables /. sol /. t -> t0},
       {{Red, Thin,
          Circle, {centers, 
           MapThread[EuclideanDistance, {centers, joints}]}]},
        Line[pts], Line[Transpose[{centers, pts}]],
        EdgeForm[Black], Red, Disk[#, 0.1] & /@ pts
     PlotRange -> {{0, wheels + 2}, {-2, 2}}
  • $\begingroup$ I think you're still missing the piston... $\endgroup$
    – rm -rf
    Jan 16, 2014 at 18:40
  • 6
    $\begingroup$ @rm-rf If I was paid... $\endgroup$
    – Michael E2
    Jan 16, 2014 at 18:43
  • $\begingroup$ Ah, yes of course! I thought you just missed it :) This should be sufficient for the OP.. $\endgroup$
    – rm -rf
    Jan 16, 2014 at 18:44
  • $\begingroup$ +1: Amazing!! That's exactly what I wanted. Thanks (I'll wait for the other answers for a while) ;-) $\endgroup$ Jan 17, 2014 at 2:18
  • $\begingroup$ @Waffle'sCrazyPeanut You're welcome! NDSolve is pretty amazing. $\endgroup$
    – Michael E2
    Jan 17, 2014 at 2:19

I did a solution with contour tracing on the distance function. It gets pretty unstable sometimes, but it's a fun question to experiment with interactivity.


DynamicModule[{p1 = {0, 2}, p2 = {1, 3}, angles = {0, 0}, distance, 
  grad, tangent}, 
 distance[a1_, a2_] := 
  Norm[{Cos@a1, Sin@a1} - (Norm[p2 - p1] {Cos@a2, Sin@a2} + p1)]; 
 grad = D[distance[a1, a2], {{a1, a2}}] /. Abs' -> Sign; 
 tangent = Cross@grad; 
    angles = 
     Mod[# /. FindRoot[distance @@ # == distance @@ angles, {t, 0}] &[
       grad t + angles + .05 tangent /. Thread[{a1, a2} -> angles]], 
      2 Pi]; Graphics[{{Dashed, Circle[], Circle[p1, Norm[p2 - p1]], 
       Locator@Dynamic@p1, Locator@Dynamic@p2}, Red, 
      Line[{{0, 0}, {Cos@#, Sin@#}, 
          p1 + Norm[p2 - p1] {Cos@#2, Sin@#2}, p1} & @@ angles]}, 
     PlotRange -> 5, Frame -> True]], 
     ContourPlot[distance[a1, a2], {a1, 0, 2 Pi}, {a2, 0, 2 Pi}, 
      ContourLabels -> True]]}]]
  • $\begingroup$ +1: That's a pretty plot. And, thanks for your experiment :P $\endgroup$ Jan 17, 2014 at 4:23

You can do this without a NDSolve by calculating the distance from the follower cranks joint to the end of the driving cranks end. Then use this distance with law of cosines to calculate the deviation angle. This is also pretty easy to implement on ANY hardware capable of doing a ArcCos and Atan2 operation (note that in c atan2 parameters are swapped).

With[{p0 = {0, 0}, p3 = {1, 3}, r1 = 1, r2 = 2, r3 = 3},
 Module[{p1, p2, s, d, a1, a2},
   p1 = r1 {Cos[a], Sin[a]};
   s = p1 - p3;
   d = Norm[s]; 
   a1 = ArcTan[s[[1]], s[[2]]];
   a2 = a1 + ArcCos[( d^2 + r2^2 - r3^2)/(2 d r2)];
   p2 = p3 + r2 {Cos[a2], Sin[a2]};
    Line[{p0, p1, p2, p3}],
    PlotRange -> {{-2, 6}, {-2, 6}}
   {a, 0, 2 Pi}]

Results in the same thing as the accepted answer:

enter image description here

This solution is much much simpler, but won't work in 3D cases. The accepted answer will, and its a more flexible model. Finally, the position of the piston, which wasn't calculated by other posters as trivial is:

pp = {r1*Cos[a] + r4 + Sin[a]^2/(2 r4), 0};

And you get the following by fiddling a bit with the graphics elements:

enter image description here

  • 2
    $\begingroup$ Nice. Of course I would say you can do this without trigonometry. :P $\endgroup$
    – Michael E2
    Apr 3, 2014 at 21:00
  • $\begingroup$ @MichaelE2 Yes, but there are quite a few known analytical solutions outlined in literature. This one is generally one of the best behaving ones for a fourbar mechanism, without making the solution unstable as mentioned by MichaelHale. Please note the piston solution is not entirely trivial trigonometric manipulation, but there exists also a pretty trivial version that uses again the law of cosines. (BTW is your intention to make a arms race?) $\endgroup$
    – joojaa
    Apr 4, 2014 at 6:44
  • $\begingroup$ Please understand my remark was in jest. When I first started programming, doing what NDSolve does would have taken prohibitively long to program. My first inclination with this problem was to approach it as you did -- it feels like the "real" solution. (My prejudice has been slow to dissolve.) On a whim, I tried NDSolve and was surprised how easy it was. A few days later I thought of the random linkage, but I didn't bother to update until your answer bumped the question to the top of the stack. $\endgroup$
    – Michael E2
    Apr 4, 2014 at 9:41
  • $\begingroup$ @MichaelE2, No offense taken, I understood that it was a reference to my previous comment should have put a :) instead of +1. I was wondering how do you do a While loop that terminates once it does not get any warnings. NOTE: I have no problems iwth NDSolve. The NDSolve is not wrong, in fact when you start to interact with forces/bending links it becomes the ONLY way. But here the system is deterministic (a textbook problem) assuming its rigid. I do quite a lot of mechanical simulation with Adams. Your prejudice is common in my students, but simulation is used because its more accurate. $\endgroup$
    – joojaa
    Apr 4, 2014 at 10:49
  • $\begingroup$ @joojaa I'm not sure what While loop you're talking about. Perhaps using WhenEvent in NDSolve? It could be used to stop integration after the driver completes one revolution, if it does. (I did not do that in my code, though.) $\endgroup$
    – Michael E2
    Apr 4, 2014 at 12:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.