### Optimal Road Problem

### by

### Daniel Chernowitz

A thesis submitted in partial fulfillment for the degree of Bachelor of Mathematics

in the

Faculty of Mathematics and Natural Sciences Department or School Name

### February 2014

### Abstract

Faculty of Mathematics and Natural Sciences Department or School Name

Bachelor of Mathematics

by Daniel Chernowitz

We consider the problem of paving a path through a mountainous terrain. An optimal path is one that balances the requirement of a route to be both as short as possible and as flat as possible. A new form of Riemannian metric, ˜S, is proposed to weigh the length of a path, such as to penalize vertical transport more than horizontal, reflecting the cost and difficulty of traffic to change altitude. A geodesic γ of this metric is then a candidate for an optimal path. A number of numerical solution techniques to calculate geodesics are explored, most notably a functional iteration using Christoffel’s symbols of the geodesic equation, and geodesic evolution to minimize geodesic curvature. Finally, statistical analyses are done on the resulting geodesics of this weighted metric through a simulated hilly landscape. The metric performs well in balancing horizontal versus vertical travel economically, however it does not prevent a high maximum steepness over the course of the path.

Without the long, windy, sleep-deprived night busses from Cusco to Ica, I wouldn’t have anything to improve upon. For this reason, I thank the infrastructure of Peru. Besides this, my gratitude to professor Waalkens who has agreed to be the second corrector, for the subjects you have taught me so far. Finally, I thank my preceptor, professor Van der Schaft. A number of times, a lucid hint dissolved the walls I was trying to run through, and I have found it pleasant to work under you.

ii

Abstract i

Acknowledgements ii

1 Introduction 1

1.1 Motivation and Formulation . . . 2

1.1.1 Road Parameters and Costs . . . 2

1.1.2 Research Question . . . 4

2 Theory 5 2.1 Parameterization of the Landscape . . . 5

2.2 Weighted Path Length . . . 7

2.3 Geodesics . . . 12

3 Discrete Solution Techniques 15 3.1 Initial Estimate . . . 15

3.1.1 Dijkstra’s Algorithm . . . 18

3.2 Refinement . . . 21

3.2.1 Functional Iteration . . . 21

3.2.2 Respacing . . . 23

3.2.3 Geodesic Evolution . . . 25

3.3 Program Structure . . . 28

3.3.1 Terrain . . . 28

3.3.2 Creating Geodesics . . . 29

4 Experiment 32 4.1 Setup . . . 32

4.2 Results. . . 32

4.2.1 Distribution . . . 33

4.2.2 Discussion. . . 35

4.2.3 Comparison to Euclidean metric . . . 37

A Appendix 38 A.1 Refining Geodesics . . . 38

Bibliography 44

iii

## Introduction

When travelling a winding road over mountainous terrain, the question of the efficiency arises. Should the engineers responisble have taken a different route, the road might have been shorter, or the average traffic might require less fuel to traverse it. Given certain demands on the use of the road, there must be an optimal road to take. We consider the case of a path confined to the surface of the terrain, thus disregarding bridges and tunnels, and we only require the road to pass through a source and target point on the map. Besides these conditions, the course is free.

Figure 1.1: A winding road in Peru. This segment does not need to pass any particular landmarks while navigating the terrain, so it’s course could conceivably be rerouted to

a shape that is more fuel efficient. It is subject to optimalisation. [1]

1

### 1.1 Motivation and Formulation

Automobiles are designed to travel horizontally. Sometimes vertical transport is in- evitable, for instance if the journey’s destination is higher than the starting point. Be- sides this, it is desirable to have as flat a course as possible, for displacing a mass against gravity constitutes performing work, which costs fuel, or resources. Even though this work is released when travelling back down, no system is 100% efficient and there will be loss. Furthermore, when vertical displacement occurs, an automobile copes best with the least possible incline. In the absurd limit, we find it cannot drive up a wall. We pro- pose that the suitability for traffic decreases as the square of the incline. Small terrain fluctuations may be appropriate, however steep hills must be discouraged when paving a way.

In this research, we will introduce a new notion of length to measure how far points are from each other on a landscape. This weighted length penalizes vertical transport.

As far as possible, we formulate it such that horizontal displacement adds to the length normally, but vertical displacement enters quadratically. Then the road between two points that minimizes this length is optimal for our purposes.

1.1.1 Road Parameters and Costs

In an attempt to quantify the demands an engineer would empose on an optimal road, we compose a formula for the total costs of placement and use. Although in reality the party responsible for construction might not need to weigh the user’s expenses as heavily as its own, here we assume the ideal situation in which they minimize total expenditure for all involved. First we set the nomenclature of quantities related to the road.

Call the length of the road L. The euclidean distance between source and target points
is ∆_{01}. When traversing in one direction, keep track of the total cumulative distance of
vertical climb, Ct, and fall, Ft. Both are positive and the difference δ01:= (Ct−F_{t}) is the
change in altitude from the source to target points. We term the unnecessary vertical
travel to be the smaller of the two, U := min(Ct, Ft). Some of these measurements are
illustrated in figure 1.2.

Say the road costs R in unit expense per unit length to pave. For the average vehicle, the drag against the direction of motion costs D per unit length travelled. In reality, this will vary with the instantaneous speed, and more sophisticated integrals may take braking and accelerating into account, but here we assume drivers always follow the speed limit exactly. Displacing the average vehicle a unit of length vertically costs an additional H, after correcting for the amount of potential energy released when returning downhill.

Figure 1.2: Illustration of quantities L, C^{t}, Ft, ∆01& δ01on a cross section of a road.

Here, U = Ft.

Finally, the road is in use for Y unit time, and the traffic consists of Q vehicles per unit time.

Definition 1.1. The total (variable) cost S of a road is a measure of how much more resources are needed to travel from source to target points than would be in the ideal situation: a straight line.

S := (L − ∆_{01}) · (R + Y QD) + Y QU H (1.1)

The amounts Y Qδ01H and ∆01(R + Y QD) are treated as an unavoidable expenses, and the quality of the road should only be judged on the amount of excess expense. When examining a particular course, we need only know the values H, R & D to a common factor, and the total number of vehicles expected to travel the road during its lifetime, Y Q.

As rough estimates, we will use the following values.

Quantity Magnitude

R 300

D 1 · 10^{−4}

H 2 · 10^{−4}

Y Q 2 · 10^{7}

These are not substantiated here, as this is not a civil engineering research. The magni- tudes are derived loosely from approximate construction and gas prices, vehicle size and performances of current day, and the first three entries may be interpreted as costs in USD per meter. More appropriate values may be chosen for any particular study, these serve merely to exemplify the ruling later in this research.

As a final measure, an interesting statistic is the maximum steepness encountered. If it is too great, no matter the quality of the rest, this bottleneck will prohibit some traffic from traversing the road. This may be for vehicles lacking the traction or power to ascend or descend safely. The maximum steepness, max ∠, is defined as the maximum ratio of vertical to horizontal displacement.

1.1.2 Research Question

The research question is threefold, we state it as follows:

• Geometrically, what mathematical conditions should be satisfied by an optimal road?

• What solution techniques may be used to calculate the form of such a road?

• Given a collection of roads constructed with different parameters, how does their quality for traffic depend on these parameters and which values minimize the ex- pression in equation 1.1?

These questions will be answered in order in chapters 2, 3 and 4, respectively.

In the following, we will use the mathematical term path when referring to the equivalent of a road in an abstract landscape.

## Theory

### 2.1 Parameterization of the Landscape

An arbitrary 3D landscape may be seen as a 2-dimensional manifold in 3-dimensional euclidean space, using cartesian coordinates. We define U :

Definition 2.1.

U (u, v) = [x, y, z] = [u, v, V (u, v)] (2.1)

U (u, v) is a 3-valued function of 2 variables. Its first and second partial derivatives are
required to exist. The boundaries of the independent plane are arbitrary, but we will
take V and thus U to be defined on an open connected set in R^{2}. This construction is
referred to as a Monge Patch. [2]

Definition 2.2. The first fundamental form of the landscape U , or Riemannian metric [3], is given by the following functions. [4]

E :=

∂U

∂u

2

=

1, 0,∂V

∂u

2

= ∂V

∂u

2

+ 1 (2.2a)

F := ∂U

∂u ·∂U

∂v =

1, 0,∂V

∂u

·

0, 1,∂V

∂v

= ∂V

∂u

∂V

∂v (2.2b)

G :=

∂U

∂v

2

=

0, 1,∂V

∂v

2

= ∂V

∂v

2

+ 1 (2.2c)

Arranged in a 2 × 2 matrix, because this is a 2-dimensional manifold, we obtain:

S :=

"

E F

F G

#

(2.3) 5

This defines a Euclidean metric, or distance function, on the (u, v)-plane, corresponding
to the augmented space between two points due to variation in height. At each point
p ∈ U , S may be used to construct an inner product between two (row) vectors l, k in
the tangent space at p, T_{p}U . [5]

hl|ki = lSk^{T} (2.4)

In particular, this naturally defines a norm on the tangent space, and thus a vector length.

|| l || =p

lSl^{T} (2.5)

Definition 2.3. We would now like to define the parameterization of a one-dimensional path γ in the (u, v)-plane.

γ(t) = (u(t), v(t))

t ∈ [a, b]

γ(a) = (x0, y0)
γ(b) = (x_{1}, y_{1})

(2.6)

We will also call (x0, y0) the source and (x1, y1) the target of the path.

This implicitely denotes a path on the 3-dimensional manifold given by the composition (U ◦ γ).

The conventional length of the path (U ◦ γ), using the metric derived in equations2.2, is the integral of the length of the velocity vector. From equation 2.5 then follows an expression for the length: [4]

L(γ) = Z b

a

E ˙u^{2}+ 2F ˙u ˙v + G ˙v^{2}^{1}_{2}
dt =

Z b a

˙
γS ˙γ^{T}^{1}_{2}

dt (2.7)

V =˙ ∂V

∂uu +˙ ∂V

∂v ˙v (2.8)

An obvious result from the Pythagorean theorem is then the equality:

L(γ) = Z b

a

˙
γS ˙γ^{T}^{1}_{2}

dt = Z b

a

˙V^{2}+ ˙u^{2}+ ˙v^{2}

^{1}

2 dt (2.9)

### 2.2 Weighted Path Length

Before we continue with optimization, we would like to weigh path length to penalize vertical travel more than horizontal travel. Then a path minimizing length, for our purposes called a geodesics, will avoid going up and down when possible. The formulation of geodesics will be made precise in the next section.

We would like to devise a new potential function W (u, v) and corresponding manifold
U_{W} = (u, v, W (u, v)) on which the path length of (U_{W} ◦ γ) is the conventional path
length, and equals the weighted length of γ in the landscape defined by V (u, v). If we
find this W , we have reduced the question to the search for geodesics on U_{W}.

Definition 2.4. The weighted length over the landscape of U_{V}, subscript for clarity, is
given by:

L_{V}(γ) =
Z b

a

( (∂V

∂u

2

+ 1) ˙u^{2}+ 2(∂V

∂u

∂V

∂v) ˙u ˙v + (∂V

∂v

2

+ 1) ˙v^{2}

^{1}_{2}

+ α( ˙V , u, v) )

dt (2.10) It is implicit that ( ˙u, ˙v) = ( ˙u(t), ˙v(t)) = ˙γ(t). There is a linear added weight term α( ˙V , u, v), which should depend on at least the speed of vertical transport ˙V , but may also involve the position on the landscape, which is assumed known.

V is a potential function of cartesian coordinates (u, v) ∈ R^{2}, so we may use the gradient
identity

∇V = ∂V

∂u,∂V

∂v

(2.11)

Then

L_{V}(γ) =
Z b

a

n

(∇V · ˙γ)^{2}+ || ˙γ||^{2}^{1}_{2}

+ α( ˙V , u, v)o

dt (2.12)

Now we impose:

L_{V}(γ) = LW(γ) =
Z b

a

γS˙ Wγ˙^{T}^{1}_{2}
dt =

Z b a

(∇W · ˙γ)^{2}+ || ˙γ||^{2}^{1}_{2}

dt (2.13)

And we would like ∇W to be an analytic function of ∇V . The corresponding α will
arise from the difference between U_{W} and U_{V}.

Theorem 2.5. UW can only take the form of a rescaling of the vertical axis in UV, i.e.

W = aV for some a ∈ R

Proof. For W to be a well defined scalar function, the vector field ∇W , its gradient, must be conservative. This is a consequence of the gradient theorem. The criterium for a vector field (X(u, v), Y (u, v)) to be conservative is equivalent to requiring its scalar curl to vanish everywhere. The scalar curl is defined as follows:

∇ × (X, Y ) = ∂X

∂v −∂Y

∂u (2.14)

Conversely, it is trivial to show that a non-conservative vector field cannot be the gradient
of a potential function. On a simply-connected region, the conservative attribute is
sufficient to guarantee the existence of such a W . [6] Thus we are looking for operations
on ∇V that preserve the conservative feature of this vector field. Call such an operation
W = (W_{1}, W2)(x1, x2), so W = ∇W and (x1, x2) = ∇V are both conservative vector
fields.

0 = ∂W2

∂u −∂W1

∂v = ∂W2

∂x_{1}

∂x1

∂u +∂W2

∂x_{2}

∂x2

∂u −∂W1

∂x_{1}

∂x1

∂v −∂W1

∂x_{2}

∂x2

∂v (2.15)

Each instance of the function W is evaluated at the point (x_{1}, x_{2}) = ∇V . Because of
the following observations

∂x1

∂u = ∂^{2}V

∂u^{2} (2.16a)

∂x2

∂v = ∂^{2}V

∂v^{2} (2.16b)

∂x_{1}

∂v = ∂x_{2}

∂u = ∂^{2}V

∂u∂v (2.16c)

Equation 2.15 must hold for all possible functions V (u, v). Furthermore, W is an ana-
lytic vector function of two variables. It’s partial derivatives cannot be derivative opera-
tors themselves. The terms ^{∂W}_{∂x}^{2}

1

∂x1

∂u and ^{∂W}_{∂x}^{2}

2

∂x2

∂v cannot cancel, although the remaining two can. This leads us to demand:

∂W_{2}

∂x1

= ∂W_{1}

∂x2

= 0 (2.17)

So

W = (W_{1}(x1), W2(x2)) (2.18)

Substituting from equation2.16c, equation2.15 reduces to

∂

∂x1

W_{1}(x1) = ∂

∂x2

W_{2}(x2) = const. (2.19)

As W1 and W2 are functions of different variables, the equality can only hold if the derivatives both equal some q ∈ R. This means the only transformation W we can make is a linear one: for instance stretching the vertical axis.

∇W = W(∇V ) = q · ∇V (2.20a)

W (u, v) = q · V (u, v) (2.20b)

For some 1 6 q ∈ R. We choose q larger than unity to penalize vertical travel, not reward it.

For W defined as such, we ask what the form of α will be. Without loss of generality,
we temporarily rescale the time variable to t^{∗} so the speed of γ(t∗) in the (u, v)-plane is
unity, thus ˙u^{2}+ ˙v^{2}= 1. Referring back to equation 2.9, the corresponding cost function
α takes the following form:

α( ˙V , q) = γS˙ _{W}γ˙^{T}^{1}_{2}

− ˙γS_{V}γ˙^{T}^{1}_{2}

= q

(q ˙V )^{2}+ 1 −
q

( ˙V )^{2}+ 1 (2.21)

In the limits ˙V = +∞ and ˙V = −∞, α( ˙V ) will go to q

(q ˙V )^{2}+ 1 −
q

( ˙V )^{2}+ 1 ≈
q

(q ˙V )^{2}−
q

( ˙V )^{2} = |a ˙V | − | ˙V | = (q − 1)| ˙V | because q > 1. The form of the penalty is
illustrated in figure 2.1.

Though searching for distance-minimizing paths on this rescaled landscape would pro- duce paths biased against vertical transport, it is not very innovative from a mathemat- ical point of view. This amounts to finding geodesics on a new euclidean surface, which has been done in many applications before.

In light of the result of theorem 2.5, we abandon the search of a landscape function on which normal path length of γ equals the weighted path length with quadratically penalized vertical transport on U . However, an alternative Riemannian metric may be proposed to calculate this weighted path length. For this metric, there may be no well defined landscape that produces it, as in equation 2.2, however, nothing fundamental prohibits us from calculating Christoffel symbols and finding geodesic curves.

−2 −1.5 −1 −0.5 0 0.5 1 1.5 2 0

1 2 3 4

V˙

α

α(q = 3)

|2 ˙V |

Figure 2.1: Cost function α as a function of ˙V . We have taken q = 3, and far from the origin α approaches 2| ˙V | asymptotically.

From the formulation of S in definition2.2, we observe that it may be represented as:

S = T (V, u, v) + I2×2=

"

V_{u}^{2} VuVv

V_{u}V_{v} V_{v}^{2}

# +

"

1 0 0 1

#

(2.22)

Using the shorthand V_{u} := ^{∂V}_{∂u} and V_{v} := ^{∂V}_{∂v}

Here T is dependent of the shape of the landscape V , observe that ˙γT ˙γ^{T} = ˙V^{2}. The
identity component I_{2×2} stems from the fact that u & v are cartesian orthogonal co-
ordinates that always contribute to the travelling speed by the Pythagorean theorem:

˙

γI ˙γ^{T} = ˙u^{2}+ ˙v^{2}.

Definition 2.6. To weigh paths against steep vertical travel, we propose the alternative
shape operator or metric: ˜S_{q} = qT^{2}+ I.

Expanding:

S˜_{q}= q

"

V_{u}^{4}+ V_{u}^{2}V_{v}^{2} V_{u}^{3}V_{v}+ V_{v}^{3}V_{u}
V_{u}^{3}Vv+ V_{v}^{3}Vu V_{v}^{4}+ V_{u}^{2}V_{v}^{2}

# +

"

1 0 0 1

#

= q||∇V ||^{2}T + I (2.23)

The motivation behind this definition is that the steepness of the landscape enters quadratically. Additionally, we may weigh this dependency through the variable a.

The scalar function ||∇V || is the absolute steepest rate of ascent/descent at each point of U . Note we will often drop the subscript q to streamline equations.

We proceed by exploring expressions of the path length of γ, using this metric.

L_{V}(γ) =
Z b

a

˙

γ ˜S ˙γ^{T}^{1}_{2}
dt =

Z b a

q||∇V ||^{2}γT ˙˙ γ^{T} + ˙γ · ˙γ^{T}^{1}_{2}
dt

= Z b

a

q||∇V ||^{2}V˙^{2}+ ˙u^{2}+ ˙v^{2}^{1}_{2}
dt

(2.24)

The second equality is due to the distributive property of matrix multiplication.

We now wish to describe the behaviour of the cost function, α( ˙V , u, v). As before, we assume the speed of γ to be unit.

α( ˙V , q, u, v) =

γ ˜˙S ˙γ^{T}

^{1}_{2}

− ˙γS ˙γ^{T}^{1}_{2}

=

q||∇V ||^{2}V˙^{2}+ 1

^{1}_{2}

− ˙V^{2}+ 1

^{1}_{2}

(2.25)

When travelling along the direction of steepest ascent/descent with unit speed, ||∇V || =

| ˙V |. Equation2.25reduces to α = q

q ˙V^{4}+ 1 −p ˙V^{2}+ 1. The shape of the cost function
is illustrated in figure2.2.

−6 −4 −2 0 2 4 6

0 10 20 30

V˙

α

q = 1 q = 2

Figure 2.2: Cost function α as a function of ˙V , for a path along the steepest ascen- t/descent. We have taken parameters q = 1 and q = 2.

We see the function is flattened around the origin but adopts a quadratic trend for higher
vertical speeds. This quadratic trend diminishes if the travel is not along the steepest
ascent. A quirk of this metric is that the minimum cost occurs at a nonzero value of ˙V ,
although the difference |α_{min}− α_{( ˙}_{V =0)}| is small by physical standards.

If the angle between the direction of steepest ascent, and the direction of unit speed travel is θ, then the magnitude of the instantaneous vertical speed is given by | ˙V | =

cos θ·||∇V ||. Inserting this into equation2.25allows us to visualize the case of movement in all directions. See figure2.3.

−^{π}_{2}

0

π 2

−4 −2 0 2 4

0 20 40

θ

||∇V ||

α

q = 2

Figure 2.3: Cost α as a function of the local maximum steepness ||∇V || and the angle θ between the direction of γ and the direction of maximum steepness, for q = 2.

### 2.3 Geodesics

As explained in chapter 1, we will be searching for paths to connect two points that in some ways minimize the path length between them. We see the path length in equation2.24is an integral over an expression involving an undetermined function. This expression is a functional, and an extreme value of the integral can be obtained through the methods of calculus of variations. In particular, the Euler-Lagrange equations solve this problem. [? ] However, in this research, we will use the theory of geodesics, which are locally distance minimizing. Formally, that means the shortest path over a manifold connecting two points lies along a geodesic curve. It can be shown that these are the same solutions as obtained by the Euler-Lagrange equations. [3]

Geodesic curves have the property that their speed, though not necessarily velocity vector, is constant and unit. Any curve can however be reparameterized to meet this requirement.

d

dt||(U ◦ γ(t))|| = 1 (2.26)

More importantly, a geodesic’s acceleration is always normal to the surface of the man- ifold. For some C ∈ R:

d^{2}

dt^{2}(U ◦ γ(t)) = CN = C^{∗} ∂V

∂u,∂V

∂v, −1

(2.27)

Normalization of surface normal N is absorbed into the second constant C^{∗}. This
condition together with that of equation2.26are tantamount to the following differential
equations:

¨

u + Γ^{u}_{uu}u˙^{2}+ 2Γ_{uv}^{u} u ˙v + Γ˙ ^{u}_{vv}˙v^{2} = 0 (2.28a)

¨

v + Γ^{v}_{uu}u˙^{2}+ 2Γ_{uv}^{v} u ˙v + Γ˙ ^{v}_{vv}˙v^{2} = 0 (2.28b)

The functions Γ^{•}_{••} are called the Christoffel symbols of the manifold. They can be
calculated by the matrix equations: [7]

"

Γ^{u}_{uu}
Γ^{v}_{uu}

#

=

"

E F

F G

#^{−1}" _{1}

2Eu

Fu−^{1}_{2}Ev

#

(2.29a)

"

Γ^{u}_{uv}
Γ^{v}_{uv}

#

=

"

E F

F G

#^{−1}"_{1}

2E_{v}

1 2Gu

#

(2.29b)

"

Γ^{u}_{vv}
Γ^{v}_{vv}

#

=

"

E F

F G

#^{−1}"

F_{v}−^{1}_{2}G_{u}

1 2Gv

#

(2.29c)

In the above, the functions E, F and G are as defined in equations 2.2. Also, we have
used the shorthand ^{∂E}_{∂u} = Eu and similar for the other components. Finally, [ • ]^{−1}
denotes the matrix inverse, given by:

"

E F

F G

#^{−1}

= 1

EG − F^{2}

"

G −F

−F E

#

(2.30) Definition 2.7. A γ(t) that solves equations2.28is a geodesic, which locally minimizes distance with respect to the Riemannian metric S.

Putting things together, we will now construct the geodesic equation for our optimal path, as measured by the weighted metric of definition2.6. Equations 2.29produce the Christoffel symbols. To this end we insert the identities:

S˜^{−1} = 1
1 + q||∇V ||^{4}

"

q||∇V ||^{2}V_{v}^{2}+ 1 −q||∇V ||^{2}VuVv

−q||∇V ||^{2}V_{u}V_{v} q||∇V ||^{2}V_{u}^{2}+ 1

#

(2.31)

Which follows from equation2.30, and:

E_{u} = 2q (2V_{u}^{3}+ V_{v}^{2}V_{u})V_{uu}+ V_{u}^{2}V_{v}V_{uv}

(2.32a)
Ev = 2q (2V_{u}^{3}+ V_{v}^{2}Vu)Vuv+ V_{u}^{2}VvVvv

(2.32b)

Fu = q (3V_{u}^{2}Vv+ V_{v}^{3})Vuu+ (3V_{v}^{2}Vu+ V_{u}^{3})Vuv

(2.32c)

Fv = q (3V_{u}^{2}Vv+ V_{v}^{3})Vuv+ (3V_{v}^{2}Vu+ V_{u}^{3})Vvv

(2.32d)

Gu = 2q (2V_{v}^{3}+ V_{u}^{2}Vv)Vuv+ V_{v}^{2}VuVuu

(2.32e)

Gv = 2q (2V_{v}^{3}+ V_{u}^{2}Vv)Vvv+ V_{v}^{2}VuVuv

(2.32f)

These are found by direct calculation. Finally, we define the matrices Γ^{u} and Γ^{v}:

Γ^{u}=

"

Γ^{u}_{uu} Γ^{u}_{uv}
Γ^{u}_{uv} Γ^{u}_{vv}

#

, Γ^{v} =

"

Γ^{v}_{uu} Γ^{v}_{uv}
Γ^{v}_{uv} Γ^{v}_{vv}

#

(2.33)

Using this substition, the geodesic equations of2.28 take the form

¨

γ = − ˙γΓ^{u}γ˙^{T}, ˙γΓ^{v}γ˙^{T}

(2.34)

Another classification of a geodesic is that it is a curve for which the geodesic curvature
κ_{g} vanishes everywhere. Classically, we may view κ_{g} as a measure of the acceleration
of γ(t) in the tangent plane to the surface. With our altered metric, the notion of the
tangent plane is not defined here, but this curvature may nonetheless be formulated in
terms of the metric and Christoffel symbols. [2] We shall find κg useful later on.

κg :=

v u u u t

det( ˜S) h

˙
γ ˜S ˙γ^{T}

i3 · ˙γ( ˙vΓ^{u}− ˙uΓ^{v}) ˙γ^{T} + ¨u ˙v − ¨v ˙u

(2.35)

## Discrete Solution Techniques

We now move on in an attempt to solve equation 2.34. As there will not be an analytic solution for the general case of an arbitrary landscape defined by V (u, v), we turn our attention to approximation techniques. If we discretize the environment, restricting γ to a succession of straight line segments passing through a finite set of points in the (u, v) plane, then it is possible to determine the optimal path under these conditions. In this chapter, we make this notion concrete.

### 3.1 Initial Estimate

In graph theory, consider a graph with weighted edges. There are a finite number of ways to travel from a starting vertex to a goal vertex, granted we forbid self-intersection.

This may also be done in our case by creating a discrete graph to describe the contin- uous landscape, to a good degree of approximation. As a first step in discretizing the environment, we define a mesh, analogous to that in [5].

Definition 3.1. Let U (u, v) be defined at least inside the rectangular domain (u, v) ∈
[u_{B}, u_{C}] × [v_{B}, v_{C}]. Then for n ∈ N>1, the mesh M_{n} is a set of (n + 1)^{2} evenly spaced
points, of this domain, positioned on a rectangular grid. We assume the boundaries
of the domain of U are sufficiently far from the target and source that expanding the
domain further does not facilitate a better path.

M_{n}= (ui, vj) =

u_{B}+(u_{C}− u_{B})

n i, v_{B}+(v_{C}− v_{B})

n j

i, j ∈ {0, 1, ..., n} (3.1)

15

This mesh represents a graph, if we treat each grid point (ui, vj) as a vertex, connected by edges to each of its direct or diagonal neighbors. The terms grid point and vertex will be synonymous in the following.

Definition 3.2. The induced graph on the domain (u, v) ∈ [u_{B}, u_{C}] × [v_{B}, v_{C}] is the
pair of sets (Mn, Bn) in which Bn consists of subsets with two elements of Mn that are
connected by an edge.

B_{n}= ((ui, vj), (u_{i+k}, v_{j+l}))

(ui, vj), (u_{i+k}, v_{j+l}) ∈ Mn

k, l{−1, 0, 1} (3.2)

An example of such a graph is shown in figure3.1.

Figure 3.1: An example of a mesh with n=6. The boundaries of the domain are
noted, as well as a point(u_{i}, v_{j}). Also, in the bottom right, the edges that originate

from any element of M6 not on the boundary are drawn.

Definition 3.3. A discrete path γ^{∗}(t) = (u^{∗}(t), v^{∗}(t)), t ∈ [a, b] is the continuous, piece-
wise differentiable function obtained by subsequently connecting neighboring elements
of Mn by straight lines, without self-intersections. Furthermore, γ^{∗}(a) = (x0, y0) &

γ^{∗}(b) = (x_{1}, y_{1}). This places a slight constraint on the possible geodesics the path can
mimic, but with a sensible choice of u_{C}, u_{B}, n etc. this should not be an obstacle. For
an illustration, see figure3.2.

A note on notation, because the time-parameterization proves unwieldy. If γ^{∗} passes
through r grid-points, then γ^{∗}_{h} ∈ M_{n} is h^{th} vertex of the path, h ∈ {1, 2, ..., r}. In
particular, γ^{∗}(a) = γ^{∗}_{1} and γ^{∗}(b) = γ^{∗}_{r}

Figure 3.2: A discrete path γ^{∗}(t) on M_{6}. Certain mesh-points on the path are
signified.

Finally, we wish for this graph to approximate traversing the landscape. We construct
a weight function on B_{n}, such that the length of the edges mimics the metric ˜S, as
formulated in equation 2.24.

Definition 3.4. The weight function ω_{q} : B_{n} → R is constructed as a best discrete
approximation to the continuous ˜S. It is given by the following formula:

ωa: ((ui, vj), (ui+k, vj+l)) 7→

q

(∆u, ∆v) ˜S(q, ¯u, ¯v) (∆u, ∆v)^{T} (3.3)

Here, ∆u and ∆v are the horizontal step length over this edge in the x and y direction, respectively. The arguments of ˜S signify that the metric is evaluated at the midpoint between (ui, vj) and (ui+k, vj+l) in the (u, v)-plane.

For a path traced over an ordered set of points on the graph, this weight function presents
a sensible discrete appoximation of the weighted path length. L^{∗}_{V} is defined to equal the
sum of the weights of the traversed edges.

3.1.1 Dijkstra’s Algorithm

Dijkstra’s algorithm is an efficient tool to find the shortest path through a grid. It
constructs a tree, so that for each vertex in M_{n} the weighted length of the shortest
connection to source vertex γ^{∗}_{1} = (x0, y0) is tabulated. When this cartography is com-
plete, say the shortest path to the target γ^{∗}_{r} = (x_{1}, y_{1}) has length L. Backtracking
from the target vertex, the direct neighbor with the shortest path to the source is ob-
viously on the path of length L. Iteratively, continuing to hop to the neighbor with
the shortest tabulated path length will eventually connect the target to the source, and
this sequence defines the optimal path γ^{∗}. [8] The process of tabulating the smallest
distance to each vertex, modified for a graph induced on a mesh, is described in the box
denoted Algorithm 1.

Result: Tabulating the shortest distance from the source vertex to each other vertex
For mesh Mn, construct two (n + 1) × (n + 1) matrices, D with all entries equal to ∞
and E , with all entries zero. D contains the tentative shortest distance to each vertex,
and E keeps track of which vertices have been visited so far. Indices begin numbering
at zero, so D_{0,0} is the top left element of D, and in general elements are called by D_{i,j}
for i, j ∈ {0, 1, ..., n}, analogously for E . It is understood that D_{i,j} and E_{i,j} correspond
to the vertex ui, uj of Mn;

Set the current element (I, J ) to correspond to the source vertex, i.e. (u_{I}, v_{J}) =
(x0, y0). Each element will only be current once. Afterwards, the element or vertex is
said to be visited ;

Additionally, set DI,J := 0, the distance from the source to itself is 0;

while The target vertex is unvisited: Ei,j 6= ∞, where (u_{i}, vj) = (x1, y1) do

For all neighbors of the current element, check whether the path coming directly from the vertex of the current element is shorter than the distance already tabulated to that neighbor: for k, l ∈ {−1, 0, 1} do

if The neighbor is on the grid, (I + k), (J + l) ∈ {0, 1, ..., n} then
δ_{d}:= min{D_{I+k,J +l}, D_{I,J}+ ωa((u_{I}, v_{J}), (u_{I+k}, v_{J +l}))};

If applicable, replace the tentative distance at the neighbor by the new,
shorter distance. D_{I+k,J +l}:= δ_{d};

end end

Mark the current element as visited. E_{I,J} := ∞;

Reset the current element (I, J ) as that corresponding to the unvisited vertex with the smallest tentative distance, such that (D + E )I,J := min

i,j {E_{i,j}+ Di,j}.;

end

Algorithm 1: Dijkstra’s algorithm on a discretized landscape

One may note that not all vertices have to be consideren by this algorithm. If all unvisited vertices have a larger tentative distance to the source than the target, a path through them cannot be shorter than one through the collection of visited vertices.

We will not go into detail proving the working of dijkstra’s algorithm, as it is a well known and frequently used tool.[8] We observe that in this implementation, there will always be multiple connecting paths amongst which to minimize weighted length. Take, for instance, one commencing with solely horizontal travel to align the path with the u-coordinate of the target, and then vertical travel to align the v-coordinate.

We do not know a priori the resulting discrete path γ^{∗}of Dijkstra’s algorithm. Thus it is
impossible to judge in the general case how well U ◦ γ^{∗} will approximate a geodesic U ◦ γ
on the same surface. We can however construct a different discrete path to stay close to
a smooth curve. The question arises whether a discrete path’s weighted length can be
comparable to that of an arbitrary, given geodesic. We claim that for a sufficiently fine
mesh, this is the case, and we will put bounds on the ratio between the path lengths.

Then, because both the geometrical as well as the combinatorial approach produce an optimal path, we know that their respective results will not lie far apart.

Theorem 3.5. In the limit n → ∞, a discrete path γ^{∗}over the mesh M_{n}sharing source
and target points with a smooth curve γ, can be made to satisfy L_{V}(γ^{∗}) 6 CLV(γ), for
a suitable constant C ∈ R.

Proof. Consider a rectangular patch P of (m + 1) × (m + 1) grid points of the mesh
that γ = (u(t), v(t)) passes through, m ∈ N>1. This patch has dimensions ^{m}_{n}([u_{C} −
u_{B}] × [v_{C} − v_{B}]). For suffiently large n, and thus sufficiently fine mesh, we may make
the approximation that inside P , γ is a straight line, and ˜S is constant. Furthermore,
without loss of generality, we confine our attention to the case where 0 6 ^{dv}_{du} = ˙v/ ˙u 6 1.

Because the edges on the grid are symmetrical in all diagonal and right-angle directions,
any scenario in which the slope of γ lies outside the interval [0, 1], in another octant,
is analogous to this one. We then construct a test path γ^{∗} as follows: On the left
edge of the patch P , take the first point of γ^{∗} on the vertex closest to where γ enters
P . Then connect this point to a vertex to the right along the u-direction, ∆v = 0 or
diagonally up to the right, ∆v 6= 0, such that within that column of gridpoints γ^{∗} and
γ are closest together. Iterate this procedure m times, until the patch is traversed. For
an illustration, see figure3.3.

Within P , because the metric ˜S is taken constant, we may rearrange the diagonal and
straight segments of γ^{∗}to be grouped together. The number of straight edges is denoted
l, the number of diagonal edges is k and as an approximation, we will assume that inside
P , γ also travels the length of k edges in the v-direction. The weighted length LV is

Figure 3.3: The smooth curve γ is locally assumed linear and runs through a small
patch P of the mesh. The construction of a test path γ^{∗} such that it remains close to

γ is illustrated.

defined in equation2.24. The ratio of the discrete to smooth lengths in the patch is then given by:

L_{V}(γ^{∗})
L_{V}(γ)

P

= Rl

0

q

(1, 0) ˜S(1, 0)^{T}dt +Rk

0

q

(1, 1) ˜S(1, 1)^{T}dt
R(l+k)

0

q

(1,_{k+l}^{k} ) ˜S(1,_{k+l}^{k} )^{T}dt

(3.4)

Here, we have parameterized both paths by their u-coordinate and have used := ^{m}_{n}(u_{C}−
u_{B}), the ∆u between grid points. The elements of ˜S may be called by ˜S =

"

s_{11} s_{12}
s12 s22

#
.
Substituting the auxiliary variable λ := (l/k) and Σ ˜S := (s_{11}+ 2s_{12}+ s_{22}), we obtain:

L_{V}(γ^{∗})
L_{V}(γ)

P

= λ√

s_{11}+p
Σ ˜S

ps_{11}(λ + 1)^{2}+ 2(λ + 1)s_{12}+ s_{22} (3.5)
In the limiting cases λ = 0 and λ = ∞, this expression equals unity, as expected. In
between it exhibits an extreme value, which defines C_{1}, itself a function of the metric.

C1( ˜S) := max

λ∈R>0

L_{V}(γ^{∗})
L_{V}(γ)

P

= 2s11(Σ ˜S) − 2(s11+ s12)p s11Σ ˜S det( ˜S)

!^{1}_{2}

(3.6)

We observe that det( ˜S) = 1 + a^{2}||∇V ||^{4}V_{v}^{2}V_{u}^{2} > 1, so this expression is well defined and
finite. We know from the triangle inequality that this ratio can never drop under unity,
so equation 3.6 must be a maximum. The metric ˜S need not exhibit a symmetry in 8
directions. Analogously to C_{1}, we define C_{i}, i ∈ {2, 3, ..., 8} to be the ratio of discrete to
smooth weighted length, for the cases where the slope of γ|_{i} lies in the other 7 octants.

Finally, to complete the proof, we define C to be the maximum value attained by C_{i}( ˜S)
over all i and the domain of U . Then for the whole smooth path γ and it’s approximating
path discrete counterpart γ^{∗}, no matter the direction:

C = max

i max

P ⊂Mn

L_{V}(γ^{∗})
L_{V}(γ|_{i})
P

⇔ L_{V}(γ^{∗}) 6 CLV(γ) (3.7)

### 3.2 Refinement

The result of Dijkstra’s Algorithm is sequence of straight paths, none of which are
necessarily geodesic segments. This forms a first estimate to a geodesic, though, and the
methods described below may be used to refine it. From here, we abandon the notion
that the points of γ^{∗}lie on the grid, as these they may be any ordered collection points in
the (u, v)-plane. The refinement methods will perturb the position of these points, and
their resulting positions may take any values in R^{2}, inside or outside [u_{C}, u_{B}] × [v_{C}, v_{B}].

We assume U (u, v) is still well defined on γ^{∗}. Also, the discrete path length L^{∗}_{V} is
generalized to weigh each line segment analogous to equation 3.3, with ˜S evaluated
halfway between the endpoints of the segment.

3.2.1 Functional Iteration

We would like the path produced to satisfy the geodesic equation. Because γ^{∗} is a
discrete collection of points in a plane, we must construct a discrete condition that
is analogous to equation 2.34. Mimicking the technique in [5], we make the following
observations. For γ^{∗}_{h}= (u_{h}, v_{h}):

˙

u_{h}= u_{h+1}− u_{h−1}

∆th+ ∆th−1

+ O(∆t^{2}) (3.8a)

˙v_{h}= v_{h+1}− v_{h−1}

∆t_{h}+ ∆t_{h−1} + O(∆t^{2}) (3.8b)

¨

u_{h}= 1

∆t_{h}+ ∆t_{h−1}

u_{h+1}− u_{h}

∆t_{h} −u_{h}− u_{h−1}

∆t_{h−1}

+ O(∆t^{2}) (3.8c)

¨

v_{h}= 1

∆t_{h}+ ∆t_{h−1}

v_{h+1}− v_{h}

∆t_{h} −v_{h}− v_{h−1}

∆t_{h−1}

+ O(∆t^{2}) (3.8d)

In the equations above, ∆t_{h} is the time-step between γ^{∗}_{h} and γ^{∗}_{h+1}. Because a geodesic
is arclength parameterized, this is also the weighted length between these points. We
assume that for each h, these are of comparable magnitude, represented by ∆t.

If one can ensure that the time-steps between all consecutive points of γ^{∗} are sufficiently
uniform and equal to ∆t, these may be divided out. We make the approximation that
higher-order terms may be truncated.

˙

u_{h}≈ u_{h+1}− u_{h−1}

2∆t := ∆u_{h}

2∆t (3.9a)

˙v_{h}≈ v_{h+1}− v_{h−1}

2∆t := ∆v_{h}

2∆t (3.9b)

¨

u_{h}≈ u_{h+1}+ u_{h−1}− 2u_{h}

∆t^{2} (3.9c)

¨

v_{h}≈ v_{h+1}+ v_{h−1}− 2v_{h}

∆t^{2} (3.9d)

We define ∆u_{h} := u_{h+1}− u_{h−1} and ∆v_{h} equivalently.

Inserting into equation 2.34, we obtain a pair of relations to satisfy for every point of
γ^{∗}, one for the u and v coordinate, respectively.

u_{h+1}+ u_{h−1}− 2u_{h}

∆t^{2} = − ∆u_{h}
2∆t,∆v_{h}

2∆t

Γ^{u} ∆u_{h}
2∆t,∆v_{h}

2∆t

T

(3.10a)
v_{h+1}+ v_{h−1}− 2v_{h}

∆t^{2} = − ∆u_{h}
2∆t,∆v_{h}

2∆t

Γ^{v} ∆u_{h}
2∆t,∆v_{h}

2∆t

T

(3.10b)

Implicitly, the matrices γ^{•} are evaluated at γ^{∗}_{h}. We eliminate the time steps and formu-
late expressions for u_{h} and v_{h}, in terms of their neighbors. For additional brevity, set
the vector ∆γ^{∗}_{h}:= (∆u_{h}, ∆v_{h}).

u_{h} = u_{h+1}+ u_{h−1}

2 +1

8(∆γ^{∗}_{h})Γ^{u}(∆γ^{∗}_{h})^{T} (3.11a)
v_{h} = v_{h+1}+ v_{h−1}

2 +1

8(∆γ^{∗}_{h})Γ^{v}(∆γ^{∗}_{h})^{T} (3.11b)

ξ^{∗}_{h} := (uh, vh) (3.11c)

N.B. this result has a different form in [5], which the author believes to be erroneous.

We employ the equations3.11 in a method of refinement.

Definition 3.6. Define the method of functional iteration on the geodesic. For γ^{∗}
consisting of r ordered points, construct the function F : R^{2r−4}→ R^{2r−4}, for these sets
are the domain of the inner, variable points of γ^{∗}.

F : (u_{h−1}, u_{h+1}, v_{h−1}, v_{h+1}) 7→ ξ^{∗}_{h} | h ∈ {2, . . . , r − 1} (3.12)

The points (u_{1}, v_{1}) = (x_{0}, y_{0}), (u_{r}, v_{r}) = (x_{1}, y_{1}) are stationary and are not explicitely
counted as arguments of F . Iteratively applying F to the discrete path will be referred
to as functional iteration.

A discretized, true geodesic ξ^{∗} with uniform time-steps forms a fixed point of this func-
tion: F (ξ^{∗}) = ξ^{∗}, as the geodesic equations 3.11are satisfied. Empirically, this method
has the effect of smoothing γ^{∗}and, importantly, reducing the magnitude of L^{∗}_{V}(γ^{∗}) with
each iteration, suggesting the path converges on a true geodesic.

3.2.2 Respacing

The stability of functional iteration above proves to rely crucially on a uniform time- step between subsequent points of the discrete path. For this reason, we develope an algorithm for spacing to interpolate a discrete path with evenly spaced points.

Definition 3.7. A discrete path γ^{∗} consisting of r points has a spacing of δ_{s}, if the
following is satisfied.

q

(u_{h+1}− u_{h}, v_{h+1}− v_{h}) ˜S(u_{h+1}− u_{h}, v_{h+1}− v_{h})^{T} = δ_{s} |∀h ∈ {1, . . . , r − 1} (3.13)

To facilitate the construction of this spacing, ˜S is evaluated at (u_{h}, v_{h}). See algorithm
2.

Result: Interpolating a discrete path γ^{∗} with spacing δ_{s}, resulting in discrete path ξ^{∗}
δs must be significantly smaller than the minimum arc length between any two

consecutive points of γ^{∗};

The first points coincide, ξ^{∗}_{1} := γ^{∗}_{1}; H := 2, the current element of ξ^{∗};
for each point of γ^{∗} : h ∈ {1, . . . , r − 1} do

Choose ξ^{∗}_{H} on the line between γ^{∗}_{h} & γ^{∗}_{h+1}, such that it is arc length δs from ξ^{∗}_{H−1}.
δs=

q

(ξ^{∗}_{H} − ξ^{∗}_{H−1}) · ˜S(ξ_{H−1}^{∗} ) · (ξ^{∗}_{H} − ξ^{∗}_{H−1})^{T}. This amounts to solving for the
positive root of a quadratic equation with positive discriminant, we elaborate on
this calculation in equation3.16;

Set the distance to the next point of γ^{∗} in the (u, v)-plane, R := ||γ^{∗}_{h+1}− γ^{∗}_{h}||;

Calculate Σ_{R}:= ||ξ^{∗}_{H} − γ^{∗}_{h}||, the distance covered by ξ^{∗} from the current point of
γ^{∗} toward the next;

while Σ_{R}< R do

H := H + 1, the next point of ξ^{∗} becomes current;

Set θ := arctan(_{u}^{v}^{h+1}^{−v}^{h}

h+1−u_{h}), correcting for quadrant such that e1:= (cos θ, sin θ) is
the unit vector pointing from γ^{∗}_{h} to γ^{∗}_{h+1}. Place ξ^{∗}_{H} in the (u,v)-plane on the
line on the line spanned by γ^{∗}_{h} & γ^{∗}_{h+1}, an arc length δs from ξ^{∗}_{H−1}:

ξ^{∗}_{H} := ξ^{∗}_{H−1}+ δ_{s}/
q

e_{1}· ˜S(ξ^{∗}_{H−1}) · e^{T}_{1};
ΣR:= ||ξ^{∗}_{H} − γ^{∗}_{h}||;

end

The last point placed, ξ^{∗}_{H}, surpassed the original path, but it’s position will be
corrected in the next iteration of the for-loop.

end

Algorithm 2: Respacing algorithm for a discrete path

It will occur that subsequent points of ξ^{∗}lie in different segments of γ^{∗}. In general, these
segments will not be parallel, so ξ^{∗} must turn the corner. The situation is illustrated in
figure3.4.

Define ∆ξ := (γ^{∗}_{h}− ξ^{∗}_{H−1}) to be the vector component ξ^{∗} must still cover in the (u, v)-
plane before γ^{∗}_{h}. Now we call ρ the euclidean distance to be covered between γ^{∗}_{h} and
ξ^{∗}_{H}), in the direction of (γ^{∗}_{h+1}. Also, e_{1}:= (cos θ, sin θ) is constructed as in algorithm 2.

The following must hold:

δ_{s}=
q

(∆ξ + ρe_{1}) · ˜S(ξ^{∗}_{H−1}) · (∆ξ + ρe_{1})^{T} (3.14)

Now, we cannot factor out ρ as is possible if all points are collinear. We employ the quadratic formula with a discriminant, Dρ.

Figure 3.4: ξ^{∗}H−1& ξ^{∗}_{H} are on opposite sides of the current point γ^{∗}_{h}

δ^{2}_{s} = (∆ξ + ρe_{1}) · ˜S · (∆ξ + ρe_{1})^{T} (3.15a)
0 = e_{1}Se˜ ^{T}_{1}ρ^{2}+ (∆ξ ˜Se^{T}_{1} + e_{1}S(∆ξ)˜ ^{T})ρ + ∆ξ ˜S(∆ξ)^{T} − δ^{2}_{s} (3.15b)
D_{ρ}= (∆ξ ˜Se^{T}_{1} + e_{1}S(∆ξ)˜ ^{T})^{2}− 4e_{1}Se˜ ^{T}_{1}(∆ξ ˜S(∆ξ)^{T} − δ_{s}^{2}) (3.15c)

ρ =

−

∆ξ ˜Se^{T}_{1} + e_{1}S(∆ξ)˜ ^{T}

+pD_{ρ}

2e_{1}Se˜ ^{T}_{1} (3.15d)

By construction, equation3.15 has real solutions, so D_{ρ}is positive. Since also distance
is constructed to be nonnegative, we select the positive root in the final line. Note that
any metric is positive definite [3], thus the denominator is positive. With the magnitude
of ρ known, we simply place ξ^{∗}_{H} displaced by ρ along the line towards γ^{∗}_{h+1}.

ξ^{∗}_{H} := γ^{∗}_{h}+ ρe_{1} (3.16)

As a final remark, uniform spacing comes at the cost of perturbing the target point γ^{∗}_{r},
by at most δ_{s}. Although this may be acceptable, one must be wary not to allow the
target point to drift further by applying the algorithm multiple times unchecked.

3.2.3 Geodesic Evolution

After functional iteration, we explore a second method due to [9] for improving the dis- crete approximation of a geodesic. It involves calculating the local geodesic curvature