# How to plot a bicycle with square wheels

I want to draw a bicycle with square wheels similar to this picture, but I can't plot the trajectory along the curve. enter link description here

(*https://pastebin.com/3UbbfG6W*)

corner[x_] :=
Module[{θ = N@SawtoothWave[{0, 2 Pi}, x/(2 Pi)]},
Piecewise[{{{-1 - Cos[θ], Sin[θ]},
0 <= θ < Pi/2}, {{-Sqrt[2] Cos[θ - Pi/4],
Sqrt[2] Sin[θ - Pi/4]},
Pi/2 <= θ < Pi}, {{1 + Sin[θ - Pi],
Cos[θ - Pi]},
Pi <= θ < 3 Pi/2}, {{2, 0}, θ >=
3 Pi/2}}] + {4 Floor[x/(2 Pi)], 0}];

frame[θ_] :=
Module[{t}, t = If[θ <= 3.5 Pi, θ, 7 Pi - θ];
Show[Graphics[{RGBColor[0.8, 0.3, 0.2],
Polygon[{corner[t], corner[t + Pi/2] - {1, 0},
corner[t + Pi] - {2, 0}, corner[t + 3 Pi/2] - {3, 0}}], Black,
PointSize[Medium], Point[corner[t]]},
PlotRange -> {{-2, 5}, {-0.1, 1.5}}],
If[θ >= 0.5 Pi,
ParametricPlot[corner[x], {x, t, 0.5 Pi+ 0},
PlotStyle -> Black], Graphics[]]]];

Manipulate[frame[θ], {θ, 0, 7 Pi}]

• See Stan Wagon, Mathematica in Action, 3e, electronic version, RollingSquare[]. Jan 16 '20 at 13:31
• I haven't bought this paper book with CD at present,How can I get the RollingSquare[] code content in the DVD of this e-book? Jan 17 '20 at 1:07
• Thanks a lot,I've found a link to download this book's CD-ROM from mathematica.stackexchange.com/questions/87284/… Jan 17 '20 at 1:55
• You're welcome. BTW, I think both the current answers are wrong. For one thing the trajectory of a corner of a square should have a 90º angle at the point of contact with the catenary surface it rolls on. (The trajectory of a point is normal to the surface at the point of contact, and the surface has a 90º corner where the vertex of the square contacts the surface.) Jan 17 '20 at 2:13
• Before I clicked the book link I thought I was about to discover that "Stan Wagon" was the name for a wagon with square wheels. Jan 21 '20 at 22:31

## 6 Answers

This question is too interesting to resist, so I'll talk about how to analyze the problem.

Take a look at sketch above. It describes an arbitrary moment during the rolling. From the kinematics view, $$P$$ is the "instant center of rotation". From the energy view, the square's center of mass $$O$$ keeps its height, thus the potential of the square doesn't change, means it must be in balance. Either way we arrive to the same conclusion that $$\overline{OP}$$ is perpendicular to the trajectory of $$O$$ (the horizontal red dash-line).

Suppose the side length of the square is $$2$$, and the equation of the questioned curve is $$\boldsymbol{r}(s):=\left(x(s),y(s)\right)$$, where parameter $$s$$ is the length of $$\overline{CP}$$, which should be equal to the arc-length of $$\overset{\mmlToken{mo}{⏜}}{C'P\,}$$ due to the slipping-less rolling. It's straightforward to see tangent vector $$\frac{\mathrm{d}\boldsymbol{r}}{\mathrm{d}s}$$ at $$P$$ is parallel to $$\overline{CP}$$, so ($$\dot{F}$$ is a short-form for $$\frac{\mathrm{d}F}{\mathrm{d}s}$$ for any $$F$$)

$$\frac{-\dot{y}}{\dot{x}}=\frac{\mathrm{length}_\overset{\rightharpoonup}{CP}}{\mathrm{length}_\overset{\rightharpoonup}{OC}}=s\implies s\dot{x}+\dot{y}=0\;\text{.}$$

Additionally, because $$s$$ is an arc-length parameter, we have

$$\dot{x}^2+\dot{y}^2=1\;\text{.}$$

We setup the coordinate frame so the trajectory of $$O$$ lies on x-axis and $$C'$$ lies on y-axis. Solving the system is a one-liner

DSolve[{
s x'[s] + y'[s] == 0,
x'[s]^2 + y'[s]^2 == 1,
x[0] == 0,
y[0] == -1
}, {x, y}, s]


$$\left\{\left\{x\to-\sinh^{-1}(s), y\to\sqrt{s^2+1}-2\right\}, \left\{x\to\sinh^{-1}(s), y\to-\sqrt{s^2+1}\right\}\right\}$$

Selecting the solution with positive $$\dot{x}$$, we have

\left\{ \begin{align} x&=\sinh^{-1}(s) \\ y&=-\sqrt{s^2+1} \\ \end{align} \right.\;\text{,}

or

Block[{\$Assumptions = x \[Element] Reals},
y == -Sqrt[1 + s^2] /. Solve[x == ArcSinh[s], s] // FullSimplify
]


i.e.

$$y=-\cosh(x)$$

At last the animation:

ClearAll[catenaryGround, origcube, point, perp]
catenaryGround =
Plot[-Cosh[x], {x, -ArcSinh[1], ArcSinh[1]}, PlotRange -> All,
AspectRatio -> Automatic] // Cases[#, _Line, Infinity] & //
First;
origcube = {
{EdgeForm[GrayLevel[0.3]], FaceForm[GrayLevel[0.9]], Cuboid[{-1, -1}, {1, 1}]},
{GrayLevel[0.3], Line[{{0, 0}, {0, -1}}]}
};
point = {EdgeForm[{Hue[0., 1., 0.66], Thick}], FaceForm[GrayLevel[0.9]], Disk[{0, 0}, .04]};
perp = Line[{{1, 0}, {1, 1}, {0, 1}}];

ClearAll[cubeTF]
cubeTF[x_] := RotationTransform[ArcTan[1, -Sinh[x]]] /* TranslationTransform[{x, 0}]

ClearAll[periodLen, totalPeriod]
periodLen = 2 ArcSinh[1];
totalPeriod = 5;

DynamicModule[{period = 1, xshift, xC = -(periodLen/2), x, tf, center, contact, bottom},
DynamicWrapper[
Deploy@Graphics[{
{EdgeForm[GrayLevel[0.3]], FaceForm[GrayLevel[0.9]], Translate[FilledCurve@catenaryGround, {(# - 1) periodLen, 0} & /@ Range[totalPeriod]]}
, Dynamic@GeometricTransformation[origcube, tf]
, {Hue[0., 1., 0.66], Dashed, InfiniteLine[{0, 0}, {1, 0}]}
, {Hue[0.54, 1., 0.66], Dashed, Line[Dynamic@{center, contact}]}
, {Hue[0.54, 1., 0.66],
Dynamic@GeometricTransformation[perp, RightComposition[
ScalingTransform[1/8 {1, 1}],
RotationTransform[Pi/2 (<|-1 -> 2, 0 -> 2, 1 -> 3|>@Sign[x])],
TranslationTransform[center]
]]}
, {GrayLevel[0], Dynamic@GeometricTransformation[perp, RightComposition[
ScalingTransform[1/10 {1, 1}],
RotationTransform[Pi/2 (<|-1 -> 1, 0 -> 1, 1 -> 0|>@Sign[x])],
TranslationTransform[{0, -1}], tf
]]}
, {Black, AbsoluteThickness[4], CapForm[None], Line[Dynamic@{bottom, contact}]}
, {Black, AbsoluteThickness[4], CapForm[None],
Line@Dynamic[Function[s, {ArcSinh[s] + xshift, -Sqrt[1 + s^2]}] /@ N[Rescale[Rescale[Range[100]], {0, 1}, Sort@{0, Sinh[x]}]]]
}
, Dynamic@Translate[point, {center, contact}]
, Text[Style["O", Italic, 12], Dynamic[center], {0, -1}]
, Text[Style["P", Italic, 12], Dynamic[contact], Dynamic@{-Sign[x] 2, 0}]
, Text[Style["C", Italic, 12], Dynamic[bottom], Dynamic@{Sign[x] 2, -1}]
}
, ImageSize -> 800, PlotRange -> {{-1, 2 totalPeriod - 1} periodLen/2 + {-1, 1} Sqrt[2], {-1, 1} Sqrt[2]}, PlotRangePadding -> None
]
,
xC = -Cos[2 Clock[Pi, 10]] // Rescale[#, {-1, 1}, {-1, 2 totalPeriod - 1} periodLen/2] &
; center = {xC, 0}
; period = Round[xC/periodLen] + 1
; xshift = (period - 1) periodLen
; x = xC - xshift
; contact = {x, -Cosh[x]} + {xshift, 0}
; tf = cubeTF[x] /* TranslationTransform[{xshift, 0}]
; bottom = tf@{0, -1}
]
]


• Another correct answer! And such a beautiful animation! Jan 21 '20 at 12:34
• @MichaelE2 Thanks! Glad to hear my animation is liked! Jan 22 '20 at 15:26

Your curves show a bit strange case when central point keeps the same altitude. However, such situation can be simulated by the following way:

R = 1;
ϕ = π/4;

center = {#, 1} &;

crn = {center@# + {-R Cos@(# + ϕ), R Sin@(# + ϕ)},
center@# + {-R Cos@(# + ϕ + π/2),
R Sin@(# + ϕ + π/2)},
center@# + {-R Cos@(# + ϕ + π),
R Sin@(# + ϕ + π)},
center@# + {-R Cos@(# + ϕ + (3 π)/2),
R Sin@(# + ϕ + (3 π)/2)}} &;

trjc = Table[center@f, {f, 0, 5 π, π/10}];
trjcorn = Table[crn@f, {f, 0, 5 π, π/10}];
Manipulate[
Show[
ListPlot[{trjc, trjcorn[[All, 1]], trjcorn[[All, 2]],
trjcorn[[All, 3]], trjcorn[[All, 4]]}, AspectRatio -> 2/(5 π),
PlotStyle -> {Gray, Red, Orange, Green, Blue}, ImageSize -> 800,
Joined -> True],
Graphics[{
Black, PointSize[0.01], Point@center@f,
PointSize[0.01], Red, Point[crn[f][[1]]], Orange,
Point[crn[f][[2]]], Green, Point[crn[f][[3]]], Blue,
Point[crn[f][[4]]],
Black, Line@crn@f, Line[crn[f][[{1, 4}]]]},
PlotRange -> {{-π/2, 5 π}, All}]
]
, {f, 0, 5 π, π/10}]


Thus, the trjc and trjcorn contain points of trajectories of center and corners correspondingly

• I need to scroll forward along the curve drawn by the function Plot[-(Cosh[Mod[x, 2 Log[Sqrt[2] + 1], -Log[Sqrt[2] + 1.]]]) + Sqrt[ 2], {x, 0, 10}],that's the point. Jan 17 '20 at 5:26
• @Gowiththewind, and what's the problem with this? Just replace my ListPlot block with your Plot and enjoy fun :) Jan 17 '20 at 5:58
• Thanks a lot,I have basically solved this problem in the following answers. Jan 17 '20 at 8:57
• The problem is that the vertices of a square-wheeled bicycle do not trace cycloids as shown in this answer. What this answer shows is the trajectory of the square on a circular bicycle wheel, which is entirely different. Replacing the ListPlot with @Go's Plot produces even worse results. Jan 18 '20 at 11:32

To make it a bit shorter

n = 4 (*number of corners*)
crn = {Cos[#/n 2 Pi], Sin[#/n 2 Pi]} & /@ Range[n];
cen = {0, 1};
v[t_] = {t, 0};(*linear velocity*)
w[t_] = -0.2 t 2 Pi;(*angular velocity*)
tmax = 10; dt = 0.1;
trj = Table[(cen + v[t] + RotationMatrix[w[t]].#) & /@ crn, {t, 0, tmax, dt}];

ListAnimate[Table[Graphics[{Gray, Polygon[t],
Black, Point[Mean[t]],(*Centre*)
Table[{Hue[i/n], Line[trj[[All, i]]]}, {i, n}],
Black, Dashed, Line[Table[Mean[t], {t, trj}]]}]
,{t, trj}]]


You can play with different polygon and velocity. For an arbitrary polygon, you have to define its corners yourself in crn.

## Cycling on an arbitrary path

n = 4 (*number of corners*)
crn = {Cos[#/n 2 Pi], Sin[#/n 2 Pi]} & /@ Range[n];
cen = {0, 1};
v[t_] = {t, -(Cosh[Mod[t, 2 Log[Sqrt[2] + 1], -Log[Sqrt[2] + 1.]]]) +   Sqrt[2]};
(*parametric linear velocity*)
w[t_] = -0.2 t 2 Pi;(*angular velocity*)
tmax = 10; dt = 0.1;
trj = Table[(cen + v[t] + RotationMatrix[w[t]].#) & /@ crn, {t, 0, tmax, dt}];

ListAnimate[Table[Graphics[{Gray,Polygon[trj[[1]]], Polygon[trj[[-1]]],
Polygon[t], Black, Point[Mean[t]],(*Centre*)
Table[{Hue[i/n], Line[trj[[All, i]]]}, {i, n}], Black, Dashed,
Line[Table[Mean[t], {t, trj}]]}, PlotRange -> {{-1, 11}, {0, 3}}], {t, trj}]]


Note that I use a PlotRange here to fix the frame. Otherwise your animation might be shaky. A related post is Why is my GIF shaky?

• I need to scroll forward along the curve drawn by the function Plot[-(Cosh[Mod[x, 2 Log[Sqrt[2] + 1], -Log[Sqrt[2] + 1.]]]) + Sqrt[ 2], {x, 0, 10}] Jan 17 '20 at 5:24
• With your friendly help, I have basically solved this problem in the following answers.Thank you very much. Jan 17 '20 at 9:07

By the following proof, we can know that the horizontal ordinate of the center of the rectangle is the same as that of the tangent point.

(x + (1 - ArcLength[-Cosh[t] + Sqrt[2], {t, 0, x}]) Cos[
ArcTan[D[-(Cosh[x]) + Sqrt[2],
x]]] + (x - (1 + ArcLength[-Cosh[t] + Sqrt[2], {t, 0, x}]) Cos[
ArcTan[D[-(Cosh[x]) + Sqrt[2], x]]] -
2*Sin[ArcTan[D[-(Cosh[x]) + Sqrt[2], x]]]) // Simplify)


The result 2x that we get shows this.

     Animate[Show[
Graphics[{Rotate[
Translate[{EdgeForm[Black], LightGray,
Rectangle[{-1, Sqrt[2] - 1}, {1, Sqrt[2] + 1}], Black,
Line[{{-1, Sqrt[2] - 1}, {1, Sqrt[2] + 1}}],
Line[{{-1, Sqrt[2] + 1}, {1, Sqrt[2] - 1}}]}, {x, 0}],
Evaluate[
Simplify[
ArcTan[D[-(Cosh[
Mod[t, 2 Log[Sqrt[2] + 1], -Log[Sqrt[2] + 1.]]]) + Sqrt[
2], t]], t >= 0] /. t -> x], {x, Sqrt[2]}]}, Axes -> True,
AxesOrigin -> {0, 0}],
Plot[-(Cosh[Mod[t, 2 Log[Sqrt[2] + 1], -Log[Sqrt[2] + 1.]]]) + Sqrt[
2], {t, 0.001, x}, AxesOrigin -> {0, 0}], AxesOrigin -> {0, 0},
PlotRange -> {{0, 10}, {0, 3}}], {x, 0, 10},
AnimationRunning -> False, AnimationDirection -> ForwardBackward,
DefaultDuration -> 2]


• (+1) Finally, a correct answer! :) Jan 18 '20 at 2:53

Silvia already gave a pretty rigorous derivation (compare this with the treatment by Hall and Wagon), so let me show how to plot the desired trajectory of the rolling polygon's corner. One could certainly modify the code I gave here for this, but I will instead adapt this solution I previously wrote in OpenGL to Mathematica:

With[{n = 4}, (* number of sides *)
DynamicModule[{ang = π/n, apo, ha, pr, sa, traj, th, road},
apo = Cos[ang]; ha = ArcSinh[Tan[ang]];
pr = N[apo ha, 35]; sa = π/2 - ang; (* some extra precision needed *)

th[x_] := -sa - Gudermannian[ha Mod[x/pr, 2, -1]] + 2 ang Floor[(pr - x)/(2 pr)];
road[x_] := -apo Cosh[ha Mod[x/pr, 2, -1]];

traj = Table[{u, {u, 0} + AngleVector[th[u]]}, {u, -6 pr, 6 pr, pr/5}];

Manipulate[Show[Plot[road[x], {x, -6 pr, 6 pr}, AspectRatio -> Automatic,
Axes -> None, Background -> RGBColor[0., 0.169, 0.212],
Frame -> True, PlotRange -> {-2, 2},
PlotStyle -> RGBColor[0.149, 0.545, 0.824]],
Graphics[{{Gray, InfiniteLine[{0, 0}, {1, 0}]},
{Directive[FaceForm[],
EdgeForm[RGBColor[0.522, 0.6, 0.]]],
RegularPolygon[{t, 0}, {1, th[t]}, n]},
{Directive[AbsolutePointSize[5],
RGBColor[0.863, 0.196, 0.184]],
Point[Select[traj, First[#] <= t &][[All, -1]]]},
{Directive[AbsolutePointSize[3],
RGBColor[0.827, 0.212, 0.51]],
Point[{t, 0}]}}]],
{t, -6 pr, 6 pr}, SaveDefinitions -> True]]]


Here's how it looks for n = 6:

• Thanks for the Hall and Wagon reference! The Fourier series section is especially interesting! Jan 22 '20 at 15:28
• @Silvia, you're welcome; the paper is definitely a fascinating read. Jan 24 '20 at 2:53

This code is written by Shenlan

rotationPoint[φ_] := Line[{Cos[#] - φ, Sin[#] + 1} & /@ Range[φ, φ + 2 π, π/2]]
f[φ_] := 1 - Sqrt[2]/2 Sec[Mod[Pi/4 - φ, Pi/2, -Pi/4]]
Manipulate[
Show[{
Graphics[
{{Red, Circle[{φ, 1}]}, rotationPoint[-φ]},
PlotRange -> {{0, 2 Pi}, {0, 2}}],
Plot[f[φ1], {φ1, 0, φ}]}], {φ, 0.001, 2 Pi, 0.1},
TrackedSymbols :> {φ}
]


• This answer is wrong. (1) The vertices do not trace cycloids, which is what they do here. (2) The secant graph is approximately right, but zoom in shows it is not exactly right; the actual curve the square rolls on is a catenary. -- Your other answer is right. Jan 18 '20 at 11:40
• This is to use the first order approximation of catenary to solve this problem. The result is acceptable. Thank you for your guidance. Jan 18 '20 at 11:43
• The square's motion is not even approximately correct. Which answer to accept is your prerogative, of course, but accepting such a wrong solution does not seem helpful to others. Jan 18 '20 at 11:51
• I accept your suggestion and change my choice.Thanks a lot. Jan 18 '20 at 12:11