The Track of a Bicycle Back Tire

Steven R. Dunbar

sdunbar@math.unl.edu

University of Nebraska-Lincoln

Lincoln, NE, 68588-0323

Reinier Bosman

bosman@science.uva.nl

Institute for Theoretical Physics

Universiteit van Amsterdam

Amsterdam, The Netherlands

Sander Nooij

semnooij@science.uva.nl

Institute for Theoretical Physics

Universiteit van Amsterdam

Amsterdam, The Netherlands


Introduction

A rider on a bicycle goes down a road, the front tire steering a path

by following a prescribed motion, or just weaving back and forth. We

can see from the tire tracks after passing through a puddle that the

back tire follows a path similar to the path of the front tire. If we

know the path of the front tire precisely, what is the path of the

back tire? A related question is: If the front tire travels some

distance, how far does the back tire travel? It seems obvious

from experience that the back tire travels a shorter path than the

front tire. Bicycle folklore says that after a long trip the back

tire will show about 10% less wear than the front tire. Is it

possible to verify the folklore?

In this article we derive and solve differential equations for the

path of the back tire, knowing the path of the front tire.

For a parametric form of the front tire path the differential

equations for the path of the back tire are a pair of coupled

nonlinear differential equations.

These equations for the back tire path are a simple example

of vector Riccati equations and we

can solve them directly in some cases, with geometrical arguments,

with ``guess-and-check'', and with perturbation

methods, i.e. Taylor expansion methods, in other cases.We'll

then derive some quantitative rules for the distance the back tire

travels compared to the front tire.The formulation and solution of

the differential equations for the

back tire path given the front tire path is an example of a

``forward'' or ``direct problem''. In contrast, an ``inverse

problem'' would be: Given the paths of both of the front tire and the

back tire, determine which was the path of the front tire. This

inverse problem is explored in the well-known article ``Which way did

the bicycle go?'', [Wagon] which leads to an elementary calculus

problem relating the tangent lines of the two paths. We use this same

relation of the tangent lines in formulating the differential

equations for the direct problem. The solution of the direct problem

will also give another way to characterize the front and back tire

paths in terms of magnitude of the oscillations. This provides

another elementary way to solve the inverse problem in special cases

when the paths are sinusoidal functions.

Investigating the dynamics of bicycles has been a favorite topic in

physics for a long time, and there are many references in the physics

education journals, [Dilavore76, Francke90, Hunt89, Jones70,

Greenslade79, Kirshner80, Lowell82], particularly with reference to the

stability of a bicycle. There are fewer references to bicycles in the

mathematics literature, among them [Wagon, Yavin97], generally

dealing with the paths of the tires. This article provides an

elementary application of coupled nonlinear differential equations to a

familiar situation. The application and methods are suitable for

classes on methods of applied mathematics and differential equations.

Model equations

Assume the path (x_f(t), y_f(t)) of the contact point of the front

tire of a bicycle is given parametrically as a function of time. The

contact point is where the tire touches the road. The problem is to

determine the path of the rear tire. Let (x_b(t), y_b(t)) be the

parametric representation of the path of the contact point of the rear

tire. Let a be the wheelbase of the bike, that is, the distance

from the front contact point to the back contact point.

As the handlebars are turned the front tire swivels on the contact

point below, so a is a constant. The dimension of a is expressed

in the

unit of distance. Likewise, the dimensions of x_b(t), y_b(t),

x_f(t) and y_f(t) are in the units of distance.

Then the primary physical fact is that the tangent vector to the path

of the rear tire always points at the contact point on the path of the

front tire. Since the wheelbase is constant, we may express this as a

pair of differential equations for the path of the rear tire contact

point

x_b(t) + a x'_b(t)/(\sqrt{(x'_b)^2 + (y'_b)^2}) = x_f(t)

y_b(t) + a y'_b(t)/(\sqrt{(x'_b)^2 + (y'_b)^2}) = y_f(t).

This pair of equations can be rearranged as

(1a) x'_b(t)/(\sqrt{(x'_b)^2 + (y'_b)^2}) = ( x_f(t) - x_b(t))/a

(1b) y'_b(t)/(\sqrt{(x'_b)^2 + (y'_b)^2}) = ( y_f(t) - y_b(t))/a.

This equation says that the unit velocity vector for the back tire is

in the direction (unit vector!) of the wheelbase of the bicycle.

This is the primary fact about bicycle paths used in [Wagon].

Taking a closer look, equation

(1) is a relation about unit vectors. There is no

way to determine the magnitude of the velocity vector ( x'_b(t),

y'_b(t))^T . With this in

view set v(t)=\sqrt{(x'_b)^2 + (y'_b)^2} as the speed of the back tire

point.Rewriting the equations for the path of the back bicycle tire

(1) we get

(2a) x'_b(t) = v(t) ( x_f(t) - x_b(t) )/a

(2b) y'_b(t) = v(t) ( y_f(t) - y_b(t) )/a.

If we express the speed of the back tire point in terms of the known

front tire speed, then we will be able to express the

differential equation for the back-tire path in terms of (x_b(t),

y_b(t)) and known quantities.

Think of the back tire velocity as though the velocity vector of the

front tire is

dragging the back tire along. The back tire will have a velocity which

is the component of the front tire velocity in the direction of the

back tire's motion. The back tire speed is the magnitude of the

projection of the velocity vector of the front tire in the direction

of the unit vector of the bike (the unit vector in the direction from

the back tire to the front tire.)

This speed is then by elementary vector analysis

(3) v(t) = \begin{pmatrix}x'_f(t)\\y'_f(t)\end{pmatrix} \cdot \begin{pmatrix}\frac{x_f(t)-x_b(t)}{a}\\\frac{y_f(t)-y_b(t)}{a}\end{pmatrix}\nonumber\\

= x'_f(t)\frac{x_f(t) - x_b(t)}{a} + y'_f(t) \frac{y_f(t) - y_b(t)}{a}.

With this, the right sides of equations (2) are

expressed in terms of the given quantities ( x_f(t), y_f(t)) and

(x'_f(t), y'_f(t)) and the unknown track (x_b(t), y_b(t)) .

Expression as a Riccati Equation

In this section, we rewrite the equations

for the bicycle back tire in

an alternative form. This alternative form relates the equations

to the standard theory for matrix Riccati equations. The application

of matrix Riccati equations

to bicycle tire tracks is a very elementary and easily derived

example, other less elementary applications of Riccati equations arise

in transmission line theory, random noise theory, variational

equations, and control theory [Reid72]. Although we will

not pursue solving the equations in this form, the matrix form

suggests some other approaches to analyzing the equations

theoretically and solving them efficiently with numerical methods.

Insert expression (3) and collect terms

(4a) x'_b(t) = (1/a^2 [ x'_f(t) { x_f(t) - x_b(t)}^2 +

y'_f(t){ x_f(t) - x_b(t)}{ y_f(t) - y_b(t)}]

(4b) y'_b(t) = (1/a^2) [ x'_f(t){ x_f(t) - x_b(t)}{ y_f(t) - y_b(t)} +

y'_f(t){y_f(t)- y_b(t)^2.

Now the quadratic form of the equation is explicit.

Put the equations in terms of the quantities

(x_b(t) - x_f(t), y_b(t) - y_f(t))

x'_b(t) - x'_f(t) = (1/a^2){ x'_f(t){ x_f(t) - x_b(t)}^2 +

y'_f(t){ x_f(t) - x_b(t)}{ y_f(t) - y_b(t)}] - x'_f(t)

y'_b(t) - y'_f(t) = (1/a^2)[ x'_f(t){ x_f(t) - x_b(t)}{ y_f(t) - y_b(t)} +

y'_f(t){y_f(t) - y_b(t)}^2] - y'_f(t).

Now let w_1(t) = x_b(t) - x_f(t) and w_2(t) = y_b(t) - y_f(t)

(5a) w_1'(t) = (1/a^2)[ x'_f(t) w_1^2(t) + y'_f(t) w_1(t) w_2(t)] - x'_f(t)

(5b) w_2'(t) = (1/a^2)[ x'_f(t) w_1(t) w_2(t) + y'_f(t) w_2^2(t)] - y'_f(t).

Let W(t)=[w_1(t),w_2(t)]^T be a 2 \by 1 matrix (or vector).

Then equations (5) can be

expressed as

W'(t) - (1/a^2)[ W(t) F ^T(t) W(t) ] + F(t) = 0

where F( t) = [x'_f(t), y'_f(t)]^T.

This is the form of the Riccati matrix differential equation

as in [equation (2.1), page

11, Reid72]

W'(t) + W(t)A(t)

+ D(t)W(t) + W(t)B(t)W(t) - C(t) = 0

with A(t) = 0; D(t) = \mathbf{0}, the 2 \by 2 zero matrix;

B(t) = (-1/a^2) F ^T(t),

a 1 \by 2 matrix; C (t) = - F (t), a 2 \by 1 matrix; and

0 is the zero matrix of the appropriate dimensions. This

demonstrates the remarkable likeness of the equations for

simple things such as a bike to complicated

things like random noise theory. Even familiar physics can be

described by matrix Riccati equations.

Some simple paths for the front tire

Simple Case I: A straight path, no turns, front-tire speed $r$

> restart;

Consider when the front path is a straight line, say along the

x-axis. Our experience tells us that if the front-tire path is a

straight line, then the back tire follows behind on the same straight

line, lagging by only the wheelbase of the bicycle. We can verify

this directly from the differential equations which are easy to solve

in this case. Set

> x_f := v_f*t; y_f := 0;

x_f := v_f*t

y_f := 0

> xp_f := diff(x_f,t); yp_f := diff(y_f,t);

xp_f := v_f

yp_f := 0

> v := xp_f*(x_f - x_b(t))/a + yp_f*(y_f - y_b(t))/a;

v := v_f*(v_f*t-x_b(t))/a

Then the equation for the x-coordinate becomes

> eqn_x := diff(x_b(t),t) = v*(x_f -x_b(t))/a;

eqn_x := diff(x_b(t),t) = v_f*(v_f*t-x_b(t))^2/(a^2...

with initial condition

> init_x := x_b(0) = -a;

init_x := x_b(0) = -a

It is t hen easy to solve, and

even easier to check that the solution is:

> soln_x := dsolve( {eqn_x,init_x}, x_b(t));

soln_x := NULL

> simplify( subs( x_b(t) = -a + v_f*t, eqn_x));

v_f = v_f

Also the equation for the y-coordinate is easy

> eqn_y := diff(y_b(t) ,t) = subs(soln_x, v*(0 - y_b(t))/a);

eqn_y := diff(y_b(t),t) = -v_f*(v_f*t-x_b(t))*y_b(t...

with initial condition

> init_y := y_b(0) = 0;

init_y := y_b(0) = 0

Of course, the only solution of ahomogeneous linear equation with $0$ as initial condition is

identically $0$,

> dsolve( {eqn_y, init_y}, y_b(t));

y_b(t) = 0

The back tire follows a straight

line, lagging by only the wheelbase of the bicycle.

Geometric analysis of the steady-state case

Take the path of the front tire as a large circle. From our

experience riding a bike when making a circle, the back tire contact

point trails behind and inside the circular front tire path. In

steady state the bike body (wheelbase) makes a constant angle to the

large circle. The symmetry of the system shows that the angle is

constant. Since the bike is of constant length and makes a constant

angle, the back tire contact point creates a concentric circular path.

The situation is depicted in the figure below.

The question now becomes ``What is the radius of

the inner circle of the back tire?''.

> restart;
with(plottools): with(plots):
front_tire_path := circle([0,0], 10, color=red):
back_tire_path := circle([0,0], 8, color=blue):
radius_vector_back := line([0,0],[4*sqrt(2), 4*sqrt(2)]):
radius_vector_front := line([0,0], [4*sqrt(2)-6/sqrt(2), 4*sqrt(2)+6/sqrt(2)]):
bicycle :=line([4*sqrt(2), 4*sqrt(2)], [4*sqrt(2)-6/sqrt(2), 4*sqrt(2)+6/sqrt(2)]):
labels := textplot([[4,2,"b"], [0,5,"c"],[3.5,8.5,"a"]]):
display(front_tire_path, back_tire_path, radius_vector_back, radius_vector_front, bicycle, labels, scaling=constrained, axes= none);

Warning, the name changecoords has been redefined

[Maple Plot]

>

Set some notation: Let c be the radius of the large circle traced by

the front tire contact point. Let a be the wheelbase of the bike.

Let b be the (unknown) radius of the inner circle traced by the back

tire. Remember that the velocity vector of the back tire contact

point will point at the contact point of the front tire. But the

velocity vector of the back tire contact point is perpendicular to the

radius vector of the back tire contact point. That, is,

the back tire contact point radius vector of length b , the

bike wheelbase vector of length a in the direction of the back tire

velocity vector, and the front tire contact point radius vector of

length c form a right triangle. Therefore, a^2 + b^2 = c^2

and so b = sqrt{ c^2 - a^2} . The ratio of the length of the paths

of the small back tire circle to the large front tire circle is

sqrt{ c^2 - a^2 }/c = sqrt{ 1 - a^2/c^2}.

A limiting case when the front tire path is a large circle occurs as

the radius of the large circle goes to infinity. Then we can consider

the front path to be a straight line, say along the

$x$-axis. According to our results, the radius of the

back-tire path also goes to infinity, that is, also a straight line.

For a formulation in terms of the differential equations say that the

front tire starts at $(c,0)$.

The frequency of rotation is omega, so the period of turn is

(2\pi/omega, and the speed of the front tire is omega c,

x_f(t) = c cos(omega t), y_f(t) = c sin(omega t)

and so equation (3) becomes

v(t) = - (comega/a) sin(omega t){ c cos( omega t ) -

x_b(t)} + (c omega/a) cos( omega t) { sin(omega t) -

y_b(t)}.

After algebraic simplification, the differential equation for the

x-coordinate of the back contact point is then

(8a) d x_b(t)/dt} = (c omega/a^2) (c cos( omega t) - x_b(t))

( sin( omega t) x_b(t) - cos( omega t) y_b(t))

(8b) d y_b(t)/dt = (c omega/a^2) (c sin( omega t) - y_b(t))

( sin( omega t) x_b(t) - cos( omega t) y_b(t))

Note that although equations (8)

look rather difficult to solve, a solution by inspection is readily possible.

The previous geometric analysis suggests that the solution should be

of the form

x_b(t) = b cos( omega t - psi) and

y_b(t) = b sin( omega t - psi)

where b is the radius, and psi = arcsin(a/c) is the phase shift

indicating the back tire trails behind the front tire. This trial

solution can be inserted into the differential equation, and combined

and expanded using standard trigonometric identities to see that it

works with b = sqrt{c^2 - a^2}.

Simple Case III: A stunt circular turn with rear tire as pivot

A stunt circular turn with rear tire as pivot}

Another limiting case occurs when the front tire path is a small

circle. Let the radius of the front tire circle go to $a$, the wheel

base of the bicycle. That is, the front tire is turned at a right

angle to the bike body and moves in a circle of radius $a$. Then

according to our formula, the radius of the back tire path will go to

$0$. This is a stunt turn, a spin or pivot on the back

tire. Physically, the back tire contact point remains motionless.

For this special case the formulation in terms of the

differential equations is easier and more direct. Say that the

front tire starts at $(a,0)$, the back tire is positioned at the

origin.

The differential equation for the $x$-coordinate of the back contact point is then

> restart;

> x_f := a*cos(omega*t); y_f := a*sin(omega*t);

x_f := a*cos(omega*t)

y_f := a*sin(omega*t)

> xp_f := diff(x_f,t); yp_f := diff(y_f,t);

xp_f := -a*sin(omega*t)*omega

yp_f := a*cos(omega*t)*omega

> v := xp_f*(x_f - x_b(t))/a + yp_f*(y_f - y_b(t))/a;

v := -sin(omega*t)*omega*(a*cos(omega*t)-x_b(t))+co...

> eqn_x := diff(x_b(t),t) = v*(x_f -x_b(t))/a;

eqn_x := diff(x_b(t),t) = (-sin(omega*t)*omega*(a*c...

> eqn_x := expand(eqn_x);

eqn_x := diff(x_b(t),t) = sin(omega*t)*omega*x_b(t)...
eqn_x := diff(x_b(t),t) = sin(omega*t)*omega*x_b(t)...

> init_x := x_b(0) = 0;

init_x := x_b(0) = 0

> eqn_y := diff(y_b(t),t) = v*(y_f -y_b(t))/a;

eqn_y := diff(y_b(t),t) = (-sin(omega*t)*omega*(a*c...

> eqn_y := expand(eqn_y);

eqn_y := diff(y_b(t),t) = sin(omega*t)^2*omega*x_b(...
eqn_y := diff(y_b(t),t) = sin(omega*t)^2*omega*x_b(...

> init_y := y_b(0) = 0;

init_y := y_b(0) = 0

> dsolve( {eqn_x,eqn_y, init_x, init_y}, {x_b(t), y_b(t)});

Warning, computation interrupted

Note that although Maple is unable to solve the equation in a reasonable amount of time (I let it run for over 200 seconds of computation time), a solution by inspection is readily possible. Certainly, x_b(t) = 0 and y_b(t) = 0 satisfy the initial conditions and the differential equations. Therefore, by the existence-uniqueness theorem, this is the unique solution.

> subs( {x_b(t) = 0, y_b(t) = 0}, [eqn_x, eqn_y]);

[diff(0,t) = 0, diff(0,t) = 0]

>

When the front path is a sine curve

Set-up and numerical analysis

> restart;

Assume that the front tire track is a sine curve. Then it appears from experience that the back tire track must also be a sinusoidal curve with a phase shift and a smaller amplitude.

Let us start by introducing the path of the front tire, written parametrically

> x_f := s*t; y_f := A_f*sin(xi*s*t);

x_f := s*t

y_f := A_f*sin(xi*s*t)

That is, let s be a speed parameter that converts time to distance.

Let xi be a spatial frequency, so that 2 Pi/xi is the

horizontal distance over which the front-tire completes one

back-and-forth oscillation and 2 Pi/(xi s) is the time over which the

front-tire completes one oscillation. Let A_f be the amplitude of

the front tire oscillation. Again, let a be the wheelbase of the

bicycle, the distance of the front-tire contact point to the back-tire

contact point. The velocity of the back tire is

> xp_f := diff(x_f,t); yp_f := diff(y_f,t);

xp_f := s

yp_f := A_f*cos(xi*s*t)*xi*s

> v := xp_f*(x_f - x_b(t))/a + yp_f*(y_f - y_b(t))/a;

v := s*(s*t-x_b(t))/a+A_f*cos(xi*s*t)*xi*s*(A_f*sin...

> eqn_x := diff(x_b(t),t) = v*(x_f -x_b(t))/a;

eqn_x := diff(x_b(t),t) = (s*(s*t-x_b(t))/a+A_f*cos...

Equation (4) for the x-coordinate of the back-tire becomes with insertion of (9)

> eqn_x := expand(eqn_x);

eqn_x := diff(x_b(t),t) = s^3*t^2/(a^2)-2*s^2*x_b(t...
eqn_x := diff(x_b(t),t) = s^3*t^2/(a^2)-2*s^2*x_b(t...

> init_x := x_b(0) = -a;

init_x := x_b(0) = -a

The corresponding equation for the y-coordinate with insertion of (9) is

> eqn_y := diff(y_b(t),t) = v*(y_f -y_b(t))/a;

eqn_y := diff(y_b(t),t) = (s*(s*t-x_b(t))/a+A_f*cos...

> eqn_y := expand(eqn_y);

eqn_y := diff(y_b(t),t) = s^2*t*A_f*sin(xi*s*t)/(a^...
eqn_y := diff(y_b(t),t) = s^2*t*A_f*sin(xi*s*t)/(a^...

> init_y := y_b(0) = 0;

init_y := y_b(0) = 0

We choose $x_b(0) = -a$, $y_b(0) = 0$ as reasonable initial

conditions, but other initial conditions are possible.

These equations will certainly be difficult to solve analytically, so we will first solve them numerically. For numerical solution, we will need to supply some values for the parameters. We will take a = 1 , which means that if we are using SI units, the wheelbase is 1 meter which is reasonable (although just a little short for most adult bicycles which have a measured wheelbase of slightly more than one meter). Take the amplitude of oscillation A_f to be 0.3 meters, the speed

s to be 5 meters per second and the frequency xi to be 1.

> example_ivp := subs( {a = 1, A_f = 0.3, s = 5, xi = 1},
{eqn_x, eqn_y, init_x, init_y}
);

example_ivp := {y_b(0) = 0, x_b(0) = -1, diff(y_b(t...
example_ivp := {y_b(0) = 0, x_b(0) = -1, diff(y_b(t...
example_ivp := {y_b(0) = 0, x_b(0) = -1, diff(y_b(t...
example_ivp := {y_b(0) = 0, x_b(0) = -1, diff(y_b(t...

> example_back_path := dsolve( example_ivp,
{x_b(t), y_b(t)},
type = numeric, output = listprocedure
);

example_back_path := [t = proc (t) option `Copyrigh...

> with(plots):

Warning, the name changecoords has been redefined

> front_path_plot := plot( subs( {a = 1, A_f = 0.3, s = 5, xi = 1},
[x_f(t), y_f(t), t = 0..4]
),
color = red
):

> back_path_plot := odeplot( example_back_path, [x_b(t), y_b(t)], 0..4, numpoints = 80, color=blue):

> ftp := textplot([15,0.25,`front tire path`],align={RIGHT}):

> btp := textplot([19,0.1,`back tire path`],align={RIGHT}):

> display( {front_path_plot, back_path_plot}, labels = ["x (meters)","y (meters)"], scaling = constrained);

[Maple Plot]

>

> display( {front_path_plot, back_path_plot, ftp, btp}, labels = ["x (meters)","y (meters)"]);

[Maple Plot]

>

The path of the back tire has the general form of a sinusoid. Note that the amplitude of oscillation of the back tire is less than the amplitude of oscillation of the front tire, and the back tire is slightly phase-shifted behind the front-tire. This seems reasonable.

A good guess at an approximate solution

Based on the numerical solution, we can make a good guess at the approximate solution. Since it is a function similar to the front tire path, we can paramtetrize the path as x_b(t) = s t -a, and y_b(t) = A_b sin(xi s t - psi) where A_b is some amplitude, and psi is a phase shift.

From the differential equation (in a slightly different form) we know

{y_f(t) - y_b(t)}/{x_f(t) - x_b(t)} = d y_b(t)/d x_b(t).

Since x_f(t) = s t and by approximation x_b(t) = s t - a, the

difference is the wheelbase a. Using the chain rule and rearranging

the equation, we get

y_b(t) = y_f(t) - a dy_b(t)/dt * dt/d x_b(t).

Inserting the guessed solution and the known front tire track into

equation (12) yields

A_b sin(xi s t - psi) = A_f sin(xi s t) - a xi s A_b cos(xi s t - psi) 1/s.

Since this equation must hold for any time, choosing t = 0 we get

the phase shift

psi = arctan(xi a)

and using t = Pi/(2 s xi), we get the amplitude

A_b = \frac{A_f}{sqrt{1 + xi^2 a^2}}.

We can plot the guess along with the numerical solution to the

equations on scaled axes to compare them.

> x_b_app := s*t - a;
psi := arctan(xi*a); # Phase shift
y_b_app := (A_f/sqrt(1 + xi^2*a^2))*sin(xi*s*t - psi);

x_b_app := s*t-a

psi := arctan(xi*a)

y_b_app := -A_f*sin(-xi*s*t+arctan(xi*a))/(sqrt(1+x...

> back_path_app_plot := plot(subs({ a= 1, A_f = 0.3, s = 5, xi = 1}, [x_b_app(t), y_b_app(t), t = 0..4]), color = green):

> gp := textplot([-1,-0.15, "guessed path"], align={LEFT}):

> display( {front_path_plot, back_path_plot, back_path_app_plot, ftp, btp, gp}, labels = ["x (meters)","y (meters)"]);

[Maple Plot]

>

Except for a short transient, the guess seems to be identical with the numerical solution, indicating that over the long run, it was a very good guess, and the back tire path is indeed a sinusoid. The question is now how to justify the guess.

The difference between the numerically computed back tire path and the

guessed solution is plotted in below and except

for the transient the difference is about a centimeter,

less than 5% of the magnitude

of the back tire oscillation. The solution A_b sin(xi s t -psi) was a very

good guess, and the back tire path is nearly a sinusoid.

Is there another way to justify the guess?

> example_back_path;

[t = proc (t) option `Copyright (c) 1993 by the Uni...

> evalf( subs({ a= 1, A_f = 0.3, s = 5, xi = 1}, y_b_app));

-.2121320343*sin(-5.*t+.7853981634)

> difference_back_path_app_plot:=
odeplot( example_back_path,
[x_b(t),
y_b(t) - .2121320343*sin(5.*t - 0.7853981634)],
0..4, numpoints = 80, color=black):

> display(difference_back_path_app_plot, labels=["x (meters)", "y (meters)"]);

[Maple Plot]

>

A regular perturbation expansion solution

The amplitude of the front tire oscillation is small compared to the other physical parameters. Furthermore, we can look at the back tire path as the amplitude of the front tire path goes to zero. As the amplitude of the front tire path approaches zero, the front tire path should go to a straight line, for which we already know both from the geometric analysis and the inspection solution of the differential equations to be the same straight line. This suggests we take the amplitude of the front tire oscillation to be a small parameter in the differential equations and make a regular perturbation expansion in that parameter. We can insert the power series expansion into the differential equation and gather

like terms. This gives a sequence of linear differential equations

that we can solve. In applied mathematics this method

is called a regular perturbation expansion. Regular perturbation is

routinely used in all application areas where a nonlinear equation

needs to be solved, at least approximately.

> eqn_x_expansion := subs( {x_b(t) = x_b0(t) + A_f*x_b1(t), y_b(t) = y_b0(t) + A_f*y_b1(t)}, eqn_x);

eqn_x_expansion := diff(x_b0(t)+A_f*x_b1(t),t) = s^...
eqn_x_expansion := diff(x_b0(t)+A_f*x_b1(t),t) = s^...
eqn_x_expansion := diff(x_b0(t)+A_f*x_b1(t),t) = s^...
eqn_x_expansion := diff(x_b0(t)+A_f*x_b1(t),t) = s^...

> eqn_x_b0 := subs( A_f = 0, eqn_x_expansion);

eqn_x_b0 := diff(x_b0(t),t) = s^3*t^2/(a^2)-2*s^2*x...

> with(DEtools, odeadvisor);

[odeadvisor]

> odeadvisor( eqn_x_b0, x_b0(t));

[[_homogeneous, `class C`], [_1st_order, _with_line...

Maple recognizes that it is dealing with a Riccati equation!

Note that this is a Riccati

equation because of the term x_{b0}(t)^2. However, the right hand

side can be easily factored and the equation solved as either a

separable equation or by inspection

(or checked!) to yield

x_{b0}(t) = s\,t - a.

> factor( rhs(eqn_x_b0));

s*(-s*t+x_b0(t))^2/(a^2)

> x_b0 := s*t - a;

x_b0 := s*t-a

Likewise we can find the leading order equation for y_{b0}(t), using

the new information

> expansion_eqn_y := subs( {y_b(t) = y_b0(t) + A_f*y_b1(t), x_b(t) = s*t - a + A_f*x_b1(t)}, eqn_y);

expansion_eqn_y := diff(y_b0(t)+A_f*y_b1(t),t) = s^...
expansion_eqn_y := diff(y_b0(t)+A_f*y_b1(t),t) = s^...
expansion_eqn_y := diff(y_b0(t)+A_f*y_b1(t),t) = s^...
expansion_eqn_y := diff(y_b0(t)+A_f*y_b1(t),t) = s^...

> eqn_y_b0 := simplify(subs( A_f = 0, expansion_eqn_y));

eqn_y_b0 := diff(y_b0(t),t) = -s*y_b0(t)/a

The leading order equation for y_{b0}(t) is

d y_{b0}(t)/dt = -s y_{b0}(t)/a

with initial condition y_{b0}(0) = 0.

> y_b0 := rhs( dsolve( {eqn_y_b0, y_b0(0) = 0}, y_b0(t)));

y_b0 := 0

Because this equation is

linear and homogeneous with 0 as initial condition, the solution

must be identically zero

y_{b0}(t) = 0.

This proves the zero order perturbation solution agrees with

equations (6) and (7). To lowest order of

approximation, the motion of the back tire following a front tire

weaving back and forth with small amplitude is a straight line.

Now the back tire expansion, equating the terms with coefficient $A_f$

and inserting the now known solutions (16) and (18), gives

d x_{b1}(t)/dt = - 2 s x_{b1}(t)/a

with initial condition $x_{b1}(0) = 0$. The result is a homogeneous linear

differential equation, and is therefore easy to solve

\begin{equation}\label{eq:xb1solution}

x_{b1}(t) = 0.

\end{equation}

> eqn_x_b1 := simplify(
subs( {x_b0(t) = -a + s*t, y_b0(t) = 0},
subs( A_f = 0,
diff( eqn_x_expansion,A_f)
)
)
);

eqn_x_b1 := diff(x_b1(t),t) = -2*s*x_b1(t)/a

> x_b1 := rhs(dsolve( {eqn_x_b1, x_b1(0) = 0}, x_b1(t)));

x_b1 := 0

Now comes the interesting expansion for y_{b1}(t) by comparing terms

with coefficient A_f and using all of the previous information about

x_{b0}, y_{b0} and x_{b1}. The equation is

\d y_{b1}(t)/dt = (s/a) ( sin( xi s t) - y_{b1}(t) )

> eqn_y_b1 := subs( A_f=0, diff(expansion_eqn_y, A_f));

eqn_y_b1 := diff(y_b1(t),t) = s^2*t*sin(xi*s*t)/(a^...

> eqn_y_b1 := simplify(subs( y_b0(t) = 0, eqn_y_b1));

eqn_y_b1 := diff(y_b1(t),t) = s*(sin(xi*s*t)-y_b1(t...

> y_b1 := rhs(dsolve( {eqn_y_b1, y_b1(0) = 0}, y_b1(t)));

y_b1 := exp(-s*t/a)*xi*a/(1+xi^2*a^2)-(a*xi*cos(xi*...

> combine(y_b1);

(exp(-s*t/a)*xi*a-a*xi*cos(xi*s*t)+sin(xi*s*t))/(1+...

The solution with initial condition y_{b1}(0) = 0 is

y_{b1}(t) = [- a xi cos(xi s t) + a xi e^{( - (s t/a)} + sin(xi s t)]/{1 + xi^2 a^2}.

It is easy to see by using standard identities that this solution can be written as

y_{b1}(t) = sin(xi s t -psi)/{\sqrt{1+\xi^2 a^2}} + sin(psi) e^{-st/a}/{\sqrt{1+\xi^2 a^2}}

with \psi as in (13).

Now we can assemble the power series expansion of the solution using

the simply computed solutions x_{b0}(t), x_{b1}(t), y_{b1}(t),

and y_{b2}(t) into (14). This expansion gives the

equation of motion of the back tire ignoring quadratic and

higher orders of the front tire oscillation amplitude.

x_b(t) = s t -a

y_b(t) = A_f [ sin(xi s t -psi)/\sqrt{1+\xi^2 a^2} + sin(psi) e^{-st/a}/{\sqrt{1+\xi^2 a^2}}]

We could do some additional work to compare terms with coefficients

A_f^2, using the known x_{b0}(t), x_{b1}(t), y_{b0}(t),

and y_{b1}(t) to derive linear equations for x_{b2}(t) and

y_{b2}(t). However, before we do that let's

stop and examine graphically what we have so far.

> back_path_perturbation_plot := plot(subs({ a= 1, A_f = 0.3, s = 5, xi = 1}, [x_b0 + A_f* x_b1, y_b0 + A_f*y_b1, t = 0..4]), color = black):

> pertsol := textplot([13,0.23,`perturbation solution`],align={RIGHT}):

We can plot the perturbation solutions (22)

parametrically to get a sense of how the back path appears.

> display( {back_path_perturbation_plot, pertsol}, labels = ["x (meters)","y (meters)"] );

[Maple Plot]

>

> ftp := textplot([15,0.3,`front tire path`],align={RIGHT}):

> btp := textplot([1,0.23,`back tire path`],align={RIGHT}):

To relate the perturbation solution to the numerical one we plotted

the graphs

> display( {front_path_plot, back_path_plot, back_path_perturbation_plot, ftp, btp, pertsol},labels = ["x (meters)","y (meters)"] );

[Maple Plot]

>

Now even including the short transient, the perturbation is identical with the numerical solution, and the back tire path is indeed a sinusoid with a exponential transient.

The difference of the

numerically computed back tire path and the perturbation solution is

plotted in below. The transient difference is nearly zero,

and the difference is less than 2% of the amplitude of the back tire

oscillation. The back tire path is very nearly a sinusoid with an

exponential transient. This explains why the guess was a good approximation

except for the transient.

> example_back_path;

[t = proc (t) option `Copyright (c) 1993 by the Uni...

> num_back_path_app := evalf(subs({ a= 1, A_f = 0.3, s = 5, xi = 1}, y_b0 + A_f*y_b1));

num_back_path_app := .1500000000*exp(-5.*t)-.150000...

> difference_back_path_app_plot:=
odeplot( example_back_path,
[x_b(t),
y_b(t) - num_back_path_app],
0..4, numpoints = 80, color=black, view=-0.005..0.005):

>

> display(difference_back_path_app_plot, labels = ["x (meters)","y (meters)"]);

[Maple Plot]

>

A general front tire path

Now assume that the front tire path is given by x_f(t) = s t and

y_f(t) = A_f f(s t) where s is a speed parameter and

f( ) is a bounded continuously differentiable function

whose maximum value is scaled to be 1. Then the amplitude of the

motion of the front tire path is a small parameter A_f.

Additionally we assume that f(0) = 0 so the motion of the

front tire starts at the origin. Again we will discover an

approximation for the motion of the back tire path by means of regular

perturbation expansion.

Looking back at the work in the previous section, we see that the

expansions for x_{b0}(t), y_{b0}(t) and x_{b1}(t) don't involve

A_f or the front tire function f. Therefore, the equations and

their solutions will be the same, yielding again x_{b0}(t) = s t -

a, y_{b0}(t) = 0, and x_{b1}(t) = 0.

The equation corresponding to

(21) for the first-order term y_{b1}(t)

using all of this information simplifies nicely to

> restart;

> x_f := s*t; y_f := A_f*f(s*t);

x_f := s*t

y_f := A_f*f(s*t)

> xp_f := diff(x_f,t); yp_f := diff(y_f,t);

xp_f := s

yp_f := A_f*D(f)(s*t)*s

> v := xp_f*(x_f - x_b(t))/a + yp_f*(y_f - y_b(t))/a;

v := s*(s*t-x_b(t))/a+A_f*D(f)(s*t)*s*(A_f*f(s*t)-y...

> eqn_x := diff(x_b(t),t) = v*(x_f -x_b(t))/a;

eqn_x := diff(x_b(t),t) = (s*(s*t-x_b(t))/a+A_f*D(f...

> eqn_x := expand(eqn_x);

eqn_x := diff(x_b(t),t) = s^3*t^2/(a^2)-2*s^2*x_b(t...
eqn_x := diff(x_b(t),t) = s^3*t^2/(a^2)-2*s^2*x_b(t...

> init_x := x_b(0) = -a;

init_x := x_b(0) = -a

> eqn_y := diff(y_b(t),t) = v*(y_f -y_b(t))/a;

eqn_y := diff(y_b(t),t) = (s*(s*t-x_b(t))/a+A_f*D(f...

> eqn_y := expand(eqn_y);

eqn_y := diff(y_b(t),t) = s^2*t*A_f*f(s*t)/(a^2)-s^...
eqn_y := diff(y_b(t),t) = s^2*t*A_f*f(s*t)/(a^2)-s^...

> init_y := y_b(0) = 0;

init_y := y_b(0) = 0

> eqn_x_expansion := subs( {x_b(t) = x_b0(t) + A_f*x_b1(t), y_b(t) = y_b0(t) + A_f*y_b1(t)}, eqn_x);

eqn_x_expansion := diff(x_b0(t)+A_f*x_b1(t),t) = s^...
eqn_x_expansion := diff(x_b0(t)+A_f*x_b1(t),t) = s^...
eqn_x_expansion := diff(x_b0(t)+A_f*x_b1(t),t) = s^...

> eqn_x_b0 := subs( A_f = 0, eqn_x_expansion);

eqn_x_b0 := diff(x_b0(t),t) = s^3*t^2/(a^2)-2*s^2*x...

> factor(eqn_x_b0);

diff(x_b0(t),t) = s*(-s*t+x_b0(t))^2/(a^2)

> x_b0 := s*t - a;

x_b0 := s*t-a

> expansion_eqn_y := subs( {y_b(t) = y_b0(t) + A_f*y_b1(t), x_b(t) = s*t - a + A_f*x_b1(t)}, eqn_y);

expansion_eqn_y := diff(y_b0(t)+A_f*y_b1(t),t) = s^...
expansion_eqn_y := diff(y_b0(t)+A_f*y_b1(t),t) = s^...
expansion_eqn_y := diff(y_b0(t)+A_f*y_b1(t),t) = s^...

> eqn_y_b0 := subs( A_f = 0, expansion_eqn_y);

eqn_y_b0 := diff(y_b0(t),t) = -s^2*t*y_b0(t)/(a^2)+...

> y_b0 := rhs( dsolve( {eqn_y_b0, y_b0(0) = 0}, y_b0(t)));

y_b0 := 0

> eqn_x_b1 := subs( {x_b0(t) = -a + s*t, y_b0(t) = 0},
subs( A_f = 0,
diff( eqn_x_expansion,A_f)
)
);

eqn_x_b1 := diff(x_b1(t),t) = -2*s^2*x_b1(t)*t/(a^2...

> x_b1 := rhs(dsolve( {eqn_x_b1, x_b1(0) = 0}, x_b1(t)));

x_b1 := 0

> eqn_y_b1 := subs( A_f=0, diff(expansion_eqn_y, A_f));

eqn_y_b1 := diff(y_b1(t),t) = s^2*t*f(s*t)/(a^2)-s^...

> eqn_y_b1 := subs( y_b0(t) = 0, eqn_y_b1);

eqn_y_b1 := diff(y_b1(t),t) = s^2*t*f(s*t)/(a^2)-s^...

> expand(eqn_y_b1);

diff(y_b1(t),t) = s*f(s*t)/a-s*y_b1(t)/a

> y_b1 := rhs(dsolve( {eqn_y_b1, y_b1(0) = 0}, y_b1(t)));

y_b1 := Int(s*f(s*u)*exp(s*u/a)/a,u = 0 .. t)*exp(-...

Then it is easy to see that since f(st) is bounded by 1,

y_{b1}(t) is bounded by (s t)/a.We can get some better bounds

on the solution by doing some deeper analysis of the convolution.

Let w= (s/a)(t- u), then

y_{b1}(t) = \int_0^{(s/a)t } f( s t - aw ) e^{- w} dw.

Then let g(w) = f(st - a w) I_{[0, (s/a) t]}(w)

where I_{[0,S]}(w) is the indicator function on [0,S] which is 1

for w \in [0,S] and 0 otherwise. Note that g(w) is

continuous since as w \to (st)/a then f(s t - a w)

\to f(0) = 0 by the assumption that f( ) was C^1

and starts at the origin, f(0) = 0. Also g(w) is

bounded by the bound on f( ) which we have assumed to be

1. Then

y_{b1}(t) = \int_0^\infinity g(w) e^{-w} dw

so that the nature of the function as a Laplace transform is clearly revealed.

Then

|y_{b1}(t)| = | \int_0^\infinity g(w) e^{-w} dw |

\le \int_0^\infinity |g(w)| e^{-w} dw

\le \int_0^\infinity e^{-w} dw = 1.

This means that the amplitude of the back tire motion never exceeds

the amplitude of the front tire motion, a reasonable conclusion.

Another formulation and interation solution

Formulation and diagrams

For a front tire path given as f(x_f), a function of the

position x_f down the road, an interesting alternative differential

equation for the back tire path results. Of course, any front tire

path given by a function y = f(x_f) could be put in parametric form

and studied using the methods of the previous sections.

However, looking at the bicycle equations in this new form provides a simple

derivation of an unusual form of nonlinear delay differential

equation that is interesting in its own right. Using the method of

successive approximations is a typical way of solving nonlinear

equations, so the bicycle equation provides a nontrivial example of

successive approximations where we have an answer to compare against.

A good principle of research in applied mathematics is to solve a

problem in two different ways whenever possible and compare the results.

Fortunately for us, both solution methods yield the same results!

As before, take a to be the wheelbase of the bicycle,

and assume the front tire starts at the origin, so f(0) = 0.

Then the back tire starts at the coordinate (-a,0). Let the path of

the back tire be given as a function g(x_b) of the projection

of the back tire x_b on the road. Then we know for example that

g(-a) = 0.

At any time, the projection of the back tire x_b is an effective

distance

\sqrt{a^2 - ( f(x_f) - g(x_b) )^2}

behind the

projection x_f of the position of the front tire contact point, see

the schematic figure before. Then

x_f = x_b + \sqrt{a^2 - ( f(x_f) - g(x_b) )^2}.

This gives an implicit relationship between x_f and x_b, assuming

that we know the back tire path g(x_b).

> restart;

> with(plottools):

> y_f := (0.3)*sin(x): # equation (9)
psi := Pi/4: # equation (13) psi = arctan(xi*a) = arctan(1*1) = Pi/4
A_b := (0.3)/sqrt(2): # A_b = A_f/sqrt( 1 + xi^2*a^2)
y_b1 := A_b *( sin(x+1 - psi) + sin(psi)*exp(-(x+1))):
y_b := y_b1: # see equations (22a) and (22b)

> back_path := plot(y_b,
x = -1.0..2.5, -0.5..0.5,
color=black
):
front_path := plot(piecewise(x < 0, 0, y_f),
x = -1.0..2.5, -0.5..0.5,
color=black
):

> front_label := plots[textplot]([2.5, subs(x=2.5, y_f), "front tire path"], align ={ABOVE, RIGHT}):
back_label := plots[textplot]([2.5, subs(x=2.5, y_b), "back tire path"], align ={BELOW, RIGHT}):

> x_b := -0.5:
x_f := 1:

> bicycle := line( [ x_b, subs(x=x_b, y_b)],
[x_f, subs(x=x_f, y_f)],
linestyle=3):

vert := line( [x_f, subs(x=x_b, y_b)],
[x_f, subs(x=x_f, y_f)],
color = red):

base := line( [x_b, subs(x=x_b, y_b)],
[x_f, subs(x=x_b, y_b)],
color=red):

> effective_distance := line([ x_b, -0.08], [x_f,-0.08]):
tick_left := line( [x_b, -0.08], [x_b, -0.06]):
tick_right := line( [x_f, -0.08], [x_f, -0.06]):
effective_distance_label := plots[textplot]( [0.05,-0.10, "horizontal length"], align={RIGHT}):

> plots[display]({front_path, back_path, front_label, back_label,
bicycle, vert, base,
effective_distance, tick_left, tick_right, effective_distance_label
},
xtickmarks=[-2,5], ytickmarks = [-1,1], labels=["",""]);

>

[Maple Plot]

>

Then using the fundamental fact that the tangent vector to the back

tire path points in the direction of the bicycle, we can write

dg(x_b)/d x_b = (f(x_f) - g(x_b))/(\sqrt{a^2 - { f(x_f) - g(x_b) }^2}).

Again, remember that x_f is implicitly in terms of x_b, so that

there is really only one independent variable x_b in the

differential equation. We can rewrite this equation slightly more

transparently if we set x_f = x_b + \Delta. Then the differential

equation is

dg(x_b)/d x_b = ( f(x_b + \Delta) - g(x_b))/(\sqrt{a^2 - (f(x_b + \Delta) - g(x_b))^2}).

Although the differential equation now looks more

normal, it still hides quite bit of difficulty. Note that the right

side contains not only the unknown function g(x_b), but it

also has an advanced argument on the right side, x_b + \Delta.

Then the differential equation is of the form of a delay-differential

equation. But this also hides an important fact, since the delay

itself depends on the unknown function g(x_b) through the

implicit relationship x_f = x_b + \sqrt{a^2 - {f(x_f) - g(x_b)}^2}$.

Therefore, this differential equation has the

unknown function appearing not only nonlinearly on the right side, but

nonlinearly in the argument of the right side too!

A typical way to solve such an equation is by the method of successive

approximations. First, approximate the effective distance \Delta

between the front and rear tire contact points of the bicycle by the

wheelbase a. Generally, this will be an overestimate for the

effective distance, but it will be close if the amplitude of the paths

is not large.

The equation then becomes

d g(x_b)/d x_b = ( f( x_b + a) - g(x_b))/( \sqrt{a^2 - { f(x_b + a) - g(x_b)}^2)}.

For a given front tire path and a known wheelbase this equation is

easy to solve numerically.

The next approximation replaces the effective distance \Delta with

\sqrt{a^2 - {f(x_b+a) - g(x_b)}^2}

which in the figure is still larger than the true

effective distance, but clearly a closer approximation. The equation becomes

d g(x_b)/d x_b = ( f( x_b + \sqrt{a^2 - { f(x_b+a) - g(x_b)}^2} ) -

g(x_b)/{ \sqrt{a^2 - {f(x_b+ \sqrt{a^2 - { f(x_b+a) - g(x_b)}^2 } ) -

g(x_b)}^2}}.

Now the differential equation is clearly more complicated since

already the unknown function appears in the argument of the right hand

side as well as nonlinearly. Nevertheless, we can still solve the

differential equation numerically. Conceptually, one could use a

simple scheme such as the Euler method, since given an initial

condition such as g(-a) = 0, and knowing the given function

f(x_f), the value of the right hand side can be calculated to

give the slope of g(x_b) at x_b = -a. Then g(-a + h) can be

estimated, and the process can be repeated. In

practice, of course one uses a more sophisticated technique such as a

Runge-Kutta method, or a multi-step predictor-corrector method.

> with(student):
a := evalf(distance([x_b, subs(x=x_b, y_b)], [x_f, subs(x=x_f, y_f)])):

> appr1_bike := line([x_b, subs(x=x_b, y_b)],
[x_b + a, subs(x=x_b+a, y_f)],
color = magenta
):
appr1_vert := line([x_b + a, subs(x=x_b, y_b)],
[x_b + a, subs(x=x_b+a, y_f)],
color = magenta
):
appr1_base := line([x_b, subs(x=x_b, y_b)],
[x_b + a, subs(x=x_b, y_b)],
color = magenta
):

> first_appr_distance := line([ x_b, -0.13], [x_b + a,-0.13]):
appr_tick_left := line( [x_b, -0.13], [x_b, -0.11]):
appr_tick_right := line( [x_b + a, -0.13], [x_b + a, -0.11]):
appr_distance_label := plots[textplot]( [0.05,-0.15, "first approximation"], align={RIGHT}):

> plots[display]({front_path, back_path, front_label, back_label,
bicycle, vert, base,
appr1_bike, appr1_vert, appr1_base,
effective_distance, tick_left, tick_right, effective_distance_label,
first_appr_distance, appr_tick_left, appr_tick_right, appr_distance_label},
xtickmarks=[-2,5], ytickmarks = [-1,1], labels=["",""]);

[Maple Plot]

>

Numerical solution of the front and back path for comparison

In this subsection as an illustration we use the iteration procedure

to calculate the path of the back tire when the front tire path is a

sine curve as before. Then we can compare the solution numerically

calculated from the coupled Riccati equations for the parametric

form to the solution calculated numerically by the iteration

procedure described above.

The zero-th order iteration approximation equation is from

d g}(x_b)/d x_b = (f( x_b + a) - g(x_b))/( \sqrt{a^2 - {f(x_b + a) - g(x_b)}^2})

For the given front tire path f(x_f) = A_f sin(x_b), with

front wheel oscillation amplitude A_f = 0.3 and wheelbase a = 1

(the same parameters as before) this equation is easy to solve numerically.

Comparing the numerical solution to the zero-th order iteration

solution gives the figure below. Already, the

zero-th order iteration

solution is visually indistinguishable from the numerical solution.

The following code is taken directly from the section "When the front path is a sine curve: Set-up and numerical analysis". I repeat it here simply so that each section can stand-alone without depending on code from other previous sections. Output is suppressed in order to keep the size of the section down.

> restart;

> x_f := s*t: y_f := c*sin(f*s*t):

> xp_f := diff(x_f,t):
yp_f := diff(y_f,t):

> v := xp_f*(x_f - x_b(t))/a + yp_f*(y_f - y_b(t))/a:

> eqn_x := diff(x_b(t),t) = v*(x_f -x_b(t))/a:

> init_x := x_b(0) = -a:

> eqn_y := diff(y_b(t),t) = v*(y_f -y_b(t))/a:

> eqn_y := expand(eqn_y):

> init_y := y_b(0) = 0:

> example_ivp := subs( {a = 1, c = 0.3, s = 5, f = 1},
{eqn_x, eqn_y, init_x, init_y}
):

> example_back_path := dsolve( example_ivp,
{x_b(t), y_b(t)},
type = numeric, output = listprocedure
):

> with(plots):

Warning, the name changecoords has been redefined

> front_path_plot := plot( subs( {a = 1, c = 0.3, s = 5, f = 1},
[x_f(t), y_f(t), t = 0..4]
),
color = red
):

> back_path_plot := odeplot( example_back_path, [x_b(t), y_b(t)], 0..4, numpoints = 80, color=blue):

> display( {front_path_plot, back_path_plot}, labels = ["x","y"]);

[Maple Plot]

>

Numerical solution of zero-th order iteration approximation

> zeroth_soln := dsolve( {
diff(g(t), t) =
(0.3*sin(t + 1) - g(t))/sqrt( 1^2 - ( (0.3)*sin(t+1) - g(t))^2) , g(-1) = 0 }, g(t), type = numeric);

zeroth_soln := proc (rkf45_x) local i, rkf45_s, out...

> zeroth :=odeplot ( zeroth_soln, t = -1..20, -1..20, color = green):

> ftp := textplot([15,0.25,`front tire path`],align={RIGHT}):

> btp := textplot([19,0.1,`back tire path`],align={RIGHT}):

> zp := textplot([20,0.2, "iteration solution path"], align={RIGTH}):

> display( {front_path_plot, back_path_plot, zeroth, ftp, btp, zp}, labels = ["x (meters)","y (meters)"]);

[Maple Plot]

>

Numerical solution of first order iteration approximation

> first_soln :=
dsolve( { diff(g(t), t) =
(0.3*sin(t + sqrt(1^2 - (0.3*sin(t+1) - g(t))^2) ) - g(t))/sqrt(1^2 - (0.3*sin(t+sqrt(1^2 - (0.3*sin(t+1) - g(t))^2)) - g(t))^2), g(-1) = 0 }, g(t), type = numeric);

first_soln := proc (rkf45_x) local i, rkf45_s, outp...

> first := odeplot( first_soln, [t, g(t)], -1..20, color = green):

> display( {front_path_plot, back_path_plot, first}, labels = ["x","y"]);

[Maple Plot]

>

Numerical solution of second order iteration approximation

> second_soln :=
dsolve( { diff(g(t), t) =
(0.3*sin(t +
sqrt(1^2 - (0.3*sin(t+
sqrt(1^2 - (0.3*sin(t+1)-g(t))^2)
)
- g(t))^2)
) - g(t))/
sqrt(1^2 - (0.3*sin(t+sqrt(1^2 - (0.3*sin(t+sqrt(1^2 - (0.3*sin(t+1) - g(t))^2)) - g(t))^2)) - g(t))^2),
g(-1) = 0
}, g(t), type = numeric);

second_soln := proc (rkf45_x) local i, rkf45_s, out...

> second := odeplot( second_soln, [t, g(t)], -1..20, color = green):

> display( {front_path_plot, back_path_plot, first}, labels = ["x","y"]);

[Maple Plot]

>

Comparison of Iteration Approximations to Parametric Numerical Approximations and Perturbation Solutions

We can also make a tabular comparison of the iteration approximations

to the parametric numerical approximations and perturbation

solutions. This gives a comprehensive view of the quality of each of

the approximate solutions. We used for the numerical solution the

parametric equations. The comparison is shown in the table below.

> example_back_path;

[t = proc (t) option `Copyright (c) 1993 by the Uni...

> xbnp := subs( example_back_path, x_b(t));
ybnp := subs( example_back_path, y_b(t));

xbnp := proc (t) local rkf45_s, outpoint, r1, r2; g...

ybnp := proc (t) local rkf45_s, outpoint, r1, r2; g...

There is one slight problem, in that we will need to discover the time when the numeric parametric solution has an x-coordinate value equal to integer values 0,1,2,...,10. Then we can find the corresponding value of the y-coordiante of the numeric parametric solution at the same time. Then we will be able to compare the y-values of the numeric parametric solution, and the approximate functional representations at the integer values of the independent variable. In order to discover these times when the x-coordinate solution crosses integer values by using repeated bisection. We record the values in an array for later reference in the table.

> Delta := 0.05:
tol := 0.000000001: # 10^{-9}
for t from 0 to 10 do
low := ( t +1 )/5 - Delta:
high := (t + 1)/5 + Delta:
while (high - low > tol) do
try0 := (low + high)/2:
if xbnp(try0) < t then: # too small since xbnp being mono inc
low := try0:
else
high := try0:
fi:
od:
t:
xbnp(try0):
numersolution[t] := ybnp(try0):
od:

> yb1 := x ->
evalf( subs(
{ a= 1, A_f = 0.3, xi = 1 }, A_f*(-a*xi*cos(xi*(x+a))+a*xi*exp(-(x+a)/a)+sin(xi*(x+a)))/(1+xi^2*a^2)));

yb1 := proc (x) options operator, arrow; evalf(subs...

> t := 't';

t := 't'

> zeroth_soln := dsolve( { diff(g(t), t) =
(0.3*sin(t + 1) - g(t))/sqrt( 1^2 - ( (0.3)*sin(t+1) - g(t))^2), g(-1) = 0 }, g(t), type = numeric, output = listprocedure);

zeroth_soln := [t = proc (t) option `Copyright (c) ...

> g0 := subs( zeroth_soln, g(t));

g0 := proc (t) local rkf45_s, outpoint, r1, r2; glo...

> first_soln :=
dsolve( { diff(g(t), t) =
(0.3*sin(t + sqrt(1^2 - (0.3*sin(t+1) - g(t))^2) ) - g(t))/sqrt(1^2 - (0.3*sin(t+sqrt(1^2 - (0.3*sin(t+1) - g(t))^2)) - g(t))^2), g(-1) = 0 }, g(t), type = numeric, output = listprocedure);

first_soln := [t = proc (t) option `Copyright (c) 1...

> g1 := subs(first_soln, g(t));

g1 := proc (t) local rkf45_s, outpoint, r1, r2; glo...

> second_soln :=
dsolve( { diff(g(t), t) =
(0.3*sin(t +
sqrt(1^2 - (0.3*sin(t+
sqrt(1^2 - (0.3*sin(t+1)-g(t))^2)
)
- g(t))^2)
) - g(t))/
sqrt(1^2 - (0.3*sin(t+sqrt(1^2 - (0.3*sin(t+sqrt(1^2 - (0.3*sin(t+1) - g(t))^2)) - g(t))^2)) - g(t))^2),
g(-1) = 0
}, g(t), type = numeric, output = listprocedure);

second_soln := [t = proc (t) option `Copyright (c) ...

> g2 := subs( second_soln, g(t));

g2 := proc (t) local rkf45_s, outpoint, r1, r2; glo...

x numerical*solution Perturbation*Solution g0 g1 g2
0 .999753947941568466e-1 .1003572180 .100968736233442657 .999463712637670836e-1 .999762391153948272e-1
1 .219338630487217552 .2191169320 .219860049584271800 .219323809645024238 .219339002715751464
2 .177561083218232912 .1771349360 .177204093154762127 .177572234009800489 .177560680905217511
3 -.119087878936783904e-1 -.127264853e-1 -.152006079860884658e-1 -.117358209728733766e-1 -.119176609736207396e-1
4 -.186008851962797705 -.1853772770 -.187801701438022883 -.185926168658709284 -.186012672834666354
5 -.186026352812203727 -.1855660549 -.186419324081898124 -.186000410272272126 -.186027637282124569
6 -.154055374475795363e-1 -.1440056600e-1 -.126837879966191363e-1 -.155461316318330940e-1 -.153983032532182505e-1
7 .170768535258378595 .1702790615 .172766915584728897 .170682667664756066 .170772420522081969
8 .198897653840345023 .1985058236 .199470889180550670 .198866954631133491 .198899042465666531
9 .453143674974485697e-1 .4426437276e-1 .429645542693142396e-1 .454351076313592728e-1 .453081298341193594e-1
10 -.150996199006574239 -.1506598804 -.153348622904449422 -.150896192264712181 -.151000740871839484

>

>

Conclusions

We have shown that if the path of the front tire of a bicycle is

specified, it is possible to derive the corresponding path of the back

tire. In some geometrically simple cases, such as a large circular

front tire path, it is possible to derive the corresponding back tire

path precisely. In some other reasonable geometric cases, such as a

sinusoidal front tire path, it is not possible to state the

corresponding back tire path precisely, but it is possible to derive

an approximating expression to any desired degree of approximation.

In this paper, we have solved the approximation equations to

first order, which seems sufficient for most purposes. In fact, the

approximation techniques are easy to apply for any reasonably general front

tire path. The only limit to being able to express the solution

analytically is the ability to evaluate a convolution, or equivalently

to solve a first-order linear differential equation. Of course, in

any case, the equations for the back-tire path can be solved

numerically.

The coupled nonlinear differential

equations for the back tire path are easy to express and fairly

easy to solve when the front tire path is given parametrically.

When the front tire path is given directly as a function of the

position down the road, the differential equations become the more

challenging form of a delay-differential equation. The ``delay'' even

depends on the solution. Nevertheless, the solutions can still be

solved approximately through successive approximations. The solutions arrived at

numerically, through regular perturbation and through successive

approximations all agree, and with about the same amount of work, so

using one or another of the techniques should be determined by the

available information or purpose of the solution.

If the front tire path is a large circle, the back-tire path is a

concentric circle. The ratio of the circumferences is

sqrt{1-a^2/c^2}. Therefore, traveling in a circle, the back tire

should experience less wear than the front tire. Likewise, if the

front tire weaves back and forth in a sinusoidal fashion, with an

amplitude A_f and spatial frequency xi,

then the back tire also weaves back and forth in a

sinusoidal fashion with reduced amplitude A_f/sqrt{1 + a^2 xi^2}.

Although it is not possible to express the arc-length of a sinusoidal

function with a simple analytic expression, the proportionality of the

expressions for the functions show that the arc-length traveled by the

back tire is proportionately less than that traveled by the front tire.

Is it possible to verify the folklore that on a long bike trip the back

tire wear is less than that of the front tire? Probably not, even

though the analysis in this article supports the folklore. Too many

other variables intervene in the reality to be modeled so simply. For

example, if the back tire inflation is less than the front tire's, it

will wear more. The style of riding, including

braking, sliding, and skidding, can affect the wear too.

However, we can get two ``practical'' consequences from the solutions.

First, presented with two intertwined sinusoidal functions, known to be the

paths of the front and back tire of a bike, we can now confidently

know that the path with the larger amplitude is the front tire, and

the path with the proportionally smaller amplitude is the back tire

path. With additional inspection, knowing that the tangent vectors

from the back tire point with fixed distance to the front tire track,

we can find the way the bicycle went. Second, the solutions of the general front

tire path case shows that the amplitude of the back tire path never

exceeds the amplitude of the front tire path, that

is, in this model the bike doesn't ``fish-tail''.

Moderately complicated nonlinear differential equations can be found in even simple everyday

experiences. Better yet, it is possible to apply several different

solution techniques to the equations, yielding solutions of

varying kinds, which provide better understanding of both the everyday

experiences and the techniques.

Bibliography

@BOOK{reid72,

AUTHOR = "William T. Reid",

TITLE = "Riccati Differential Equations",

PUBLISHER = "Academic Press",

ADDRESS = "New York and London",

YEAR = "1972",

SERIES = "Mathematics in Science and Engineering",

VOLUME = "86"

}

@Article{kirshner80,

author = {Daniel Kirshner},

title = {Some nonexplanations of bicycle stability},

journal = {American Journal of Physics},

year = 1980,

volume = 48,

pages = {36-38}

}

@Article{jones70,

author = {David E. H. Jones},

title = {The stability of the bicycle},

journal = {Physics Today},

year = 1970,

volume = 23,

number = 04,

pages = {34-40}

}

@Article{greenslade79,

author = { Thomas B. Greenslade Jr.},

title = {Exponential bicycle gearing},

journal = {The Physics Teacher},

year = 1979,

volume = 17,

pages = {455-456}

}

@Article{greenslade83,

author = { Thomas B. Greenslade Jr.},

title = {More bicycle physics},

journal = { The Physics Teacher },

year = 1983,

volume = 21,

pages = {360-363}

}

@Article{hunt89,

author = {Robert G. Hunt},

title = {Bicycles in the physics lab},

journal = {The Physics Teacher},

year = 1989,

volume = 27,

pages = {160-165}

}

@Book{dilavore76,

author = { P. DiLavore},

title = {The Bicycle; Teacher's Guide to the Bicycle},

publisher = {(American Association of Physics Teachers},

year = 1976,

address = { College Park, MD}

}

@Article{yavin97,

author = {Y. Yavin},

title = {Navigation and control of the motion of a riderless bicycle by using a

simplified dynamic model},

journal = {Math. Comput. Modelling},

year = 1997,

volume = 25,

number = 11,

pages = {67--74},

note = {98e:70022}

}

@Comment Summary: "This work deals with the guidance and control of a riderless

@Comment bicycle. Given are two points $P\sb 1$ and

@Comment $P\sb 2$

@Comment in the horizontal plane and a finite time interval

@Comment $[0,t\sb {\rm f}]$.

@Comment Denote by $(x\sb 1,y\sb

@Comment 1,z\sb 1)$ the coordinates of the

@Comment center of the bicycle's rear wheel. Based on a simplified dynamical

@Comment model of the bicycle, and by using the concept of path

@Comment controllability, control laws are derived for the bicycle's pedalling

@Comment moment and directional moment such that $(x\sb

@Comment 1,y\sb 1)$

@Comment will move from $P\sb 1$ to

@Comment $P\sb 2$ during the time interval

@Comment $[0,t\sb {\rm f}]$."

@Article{francke90,

author = {G. Francke and W. Suhr and F. Riess},

title = {An advanced model of bicycle dynamics},

journal = {European J. Phys.},

year = 1990,

volume = 11,

number = 2,

pages = {116-121},

note = {90m:70009}

}

@Comment Summary: "A theoretical model of a moving bicycle is presented for

@Comment arbitrary bicycle geometries at finite angles. The nonlinear

@Comment equations of motion are derived and solved with the help of a

@Comment computer. The solutions are tested for energy conservation, and

@Comment examined with respect to inherent stability. For common bicycles,

@Comment velocity and lean angle ranges of self-stable motion are

@Comment predicted."

@Article{lowell82,

author = {J. Lowell and H.D. McKell},

title = {The stability of bicycles},

journal = {Amer. J. Phys.},

year = 1982,

volume = 50,

number = 12,

pages = {1106--1112},

note = {84b:70011}

}

@Comment Authors' summary: "We review previous ideas on why a bicycle is stable

@Comment and analyze a somewhat simplified model. We show

@Comment that a bicycle is almost, though not quite, self-stable; the most

@Comment important parameter governing the stability is the `castor' of the

@Comment front wheel, as was suggested by D. E. H. Jones [Phys. Today 23

@Comment (1970), no. 4, 34 - 40]. The almost stable behavior explains

@Comment why a bicycle is so easy to ride at speed."

@Book{wagon,

author = {Joseph D. E. Konhauser and Dan Velleman and Stan Wagon},

title = {Which way did the bicycle go?},

publisher = {MAA},

year = {1996}

}