Mathematical Modeling, Proof of Concept

i

Mathematical Modeling

Gerhard Dangelmayr and Michael Kirby

Department of Mathematics, Colorado State University, Fort Collins, Colorado 80524

E-mail address: gerhard@math.colostate.edu, kirby@math.university.edu

This project is supported in part by NSF Grant DUE-#0126650. \begin{matrix}c \\ ◯ \end{matrix}2002.

Contents

Preface vii

Chapter 1. Mathematical Modeling 1 1.1. Examples of Modeling 1 1.2. The Modeling Process 3 1.3. Purpose of this Course 4

Chapter 2. Qualitative Modeling with Functions 7 2.1. Modeling Species Propagation 7 2.2. Supply and Demand 8 2.3. Modeling with Proportion and Scale 11 2.4. Dimensional Analysis 18

Chapter 3. Data Fitting 27 3.1. Linear Least Squares 27 3.2. Splines 33

Chapter 4. Difference Equations 39 4.1. Analytical Solutions 41 4.2. Population Modeling 51 4.3. Empirical Modeling 54

Chapter 5. Linear Programming 65 5.1. Examples of Linear Programs 66 5.2. Geometric Solution of a 2D Linear Program 67 5.3. Sensitivity Analysis 68 5.4. Data Fitting and the Uniform Approximation 70

CHAPTER 5 Nonlinear Programming

A nonlinear programming problem has the form

min f(x)

subject to the inequality constraints

g_{i}(x) ≤ 0

for i = 1, …, p and the equality constraints

h_{i}(x) = 0

for i = 1, …, q.

Example 5.1. Given a fixed area of cardboard A construct a box of maximum volume. The nonlinear program for this is

min xyz

subject to

2xy + 2xz + 2yz = A

Example 5.2. Consider the problem of determining locations for two new high schools in a set of P subdivisions N_{j}. Let w_{1j} be the number of students going to school A and w_{2j} be the number of students going to school B from subdivision N_{j}. Assume that the student capacity of school A is c_{1} and the capacity of school B is c_{2} and that the total number of students in each subdivision is r_{j}. We would like to minimize the total distance traveled by all the students given that they may attend either school A or B. It is possible to construct a nonlinear program to determine the locations (a, b) and (c, d) of high schools A and B, respectively assuming the location of each subdivision N_{i} is model as a single point denoted (x_{i}, y_{i}).

min ∑_{j=1}^{P} w_{1j}((a − x_{j})^{2} + (b − y_{j})^{2})^{\frac{1}{2}} + w_{2j}((c − x_{j})^{2} + (d − y_{j})^{2})^{\frac{1}{2}}

subject to the constraints

∑_{j} w_{ij} ≤ c_{j}

w_{1j} + w_{2j} = r_{j}

for j = 1, …, P.

Example 5.3. Neural networks have provided a new tool for approximating functions where the functional form is unknown. The approximation takes on the form

f(x) = b_{j}σ(a_{j}x − α_{j}) − β

and the corresponding sum of squares error term is

E(a_{j}, b_{j}, α_{j}, β) = ∑_{i}(y_{i} − f(x_{i}))^{2}

The problem of minimizing the error function is, in this instance, an unconstrained optimization problem. An efficient means for computing the gradient of E is known as the backpropogation algorithm.

5.1. Unconstrained Optimization in One Dimension

Here we begin by considering a significantly simplified (but nonetheless important) nonlinear programming problem, i.e., the domain and range of the function to be minimized are one-dimensional and there are no constraints. A necessary condition for a minimum of a function was developed in calculus and is simply

f′(x) = 0

Note that higher derivative tests can determine whether the function is a max or a min, or the value f(x + δ) may be compared to f(x).

Note that if we let

g(x) = f′(x)

then we may convert the problem of finding a minimum or maximum of a function to the problem of finding a zero.

5.1.1. Bisection Algorithm. Let x^{∗} be a root, or zero, of g(x), i.e., g(x^{∗}) = 0. If an initial bracket [a, b] is known such that x^{∗} ∈ [a, b], then a simple and robust approach to determining the root is to bisect this interval into two intervals [a, c] and [c, b] where c is the midpoint, i.e.,

c = \frac{a + b}{2}

If

g(a)g(c) < 0

then we conclude

x^{∗} ∈ [a, c]

while if

g(b)g(c) < 0

then we conclude

x^{∗} ∈ [b, c]

This process may now be iterated such that the size of the bracket (as well as the actual error of the estimate) is being divided by 2 every iteration.

5.1.2. Newton’s Method. Note that in the bisection method the actual value of the function g(x) was only being used to determine the correct bracket for the root. Root finding is accelerated considerably by using this function information more effectively.

For example, image we were seeking the root of a function that was a straight line, i.e., g(x) = ax + b and our initial guess for the root was x_{0}. If we extend this straight line from the point x_{0} it is easy to determine where it crosses the axis, i.e.,

ax_{1} + b = 0

so x_{1} = −b/a. Of course, if the function were truly linear then no first guess would be required. So now consider the case that g(x) is nonlinear but may be approximated locally about the point x_{0} by a line. Then the point of intersection of this line with the x-axis is an estimate, or second guess, for the root x^{∗}. The linear approximation comes from Taylor’s theorem, i.e.,

g(x) = g(x_{0}) + g′(x_{0})(x − x_{0}) + \frac{1}{2}g′′(x_{0})(x − x_{0})^{2} + …

So the linear approximation to g(x) about the point x_{0} can be written

l(x) = g(x_{0}) + g′(x_{0})(x − x_{0})

If we take x_{1} to be the root of the linear approximation we have

l(x_{1}) = 0 = g(x_{0}) + g′(x_{0})(x_{1} − x_{0})

Solving for x_{1} gives

x_{1} = x_{0} − \frac{g(x_{0})}{g′(x_{0})}

or at the nth iteration

x_{n+1} = x_{n} − \frac{g(x_{n})}{g′(x_{n})}

The iteration above is for determining a zero of a function g(x). To determine a maximum or minimum value of a function f we employ condition that f′(x) = 0. Now the iteration is modified as as

x_{n+1} = x_{n} − \frac{f′(x_{n})}{f′′(x_{n})}

5.2. Unconstrained Optimization in Higher Dimensions

Now we consider the problem of minimizing (or maximizing) a scalar function of many variables, i.e., defined on a vector field. We consider the extension of Newton’s method presented in the previous section as well as a classical approach known as steepest descent.

5.2.1. Taylor Series in Higher Dimensions. Before we extend the search for extrema to higher dimensions we present Taylor series for functions of more than one domain variable. To begin, the Taylor series for a function of two variables is given by

g(x, y) =g(x^{(0)}, y^{(0)}) + \frac{∂g}{∂x}(x − x^{(0)}) + \frac{∂g}{∂y} (y − y^{(0)})

+ \frac{∂^{2}g}{∂x^{2}} \frac{(x − x^{(0)})^{2}}{2} + \frac{∂^{2}g}{∂y^{2}} \frac{(y − y^{(0)})^{2}}{2} + \frac{∂^{2}g}{∂x∂y} (x − x^{(0)})(y − y^{(0)})

+ higher order terms

In n variables x = (x_{1}, …, x_{n}) the Taylor series expansion becomes

g(x) = g(x^{(0)}) + ∇g(x^{(0)})(x − x^{(0)}) + \frac{1}{2}(x − x^{(0)})^{T} Hg(x^{(0)})(y − y^{(0)}) + ⋯

where the Hessian matrix is defined as

(Hg(x^{(0)}))_{ij} = \frac{∂^{2}g}{∂x_{i}∂x_{j}} |_{x=x^{(0)}}

and the gradient is written as a row vector, i.e.,

(∇g(x^{(0)}))_{i} = \frac{∂g}{∂x_{i}} |_{x=x^{(0)}}

5.2.2. Roots of a Nonlinear System. We saw that Newton’s method could be used to develop an iteration for determining the zeros of a scalar function. We can extend those ideas for determining roots of the nonlinear system

g_{1}(x_{1}, …, x_{n}) = 0
g_{2}(x_{1}, …, x_{n}) = 0
⋮
g_{n}(x_{1}, …, x_{n}) = 0

or, more compactly,

g(x) = 0.

Now we apply Taylor’s theorem to each component g_{i}(x_{1}, …, x_{n}) individually, i.e., retaining only the first order terms we have the linear approximation to g_{i} about the point x^{(0)} as

l_{i}(x) = g_{i}(x^{(0)}) + ∇g_{i}(x^{(0)})(x − x^{(0)})

for i = 1, …, n. We can write these components together as a vector equation

l(x) = g(x^{(0)}) + ∇g(x^{(0)})(x − x^{(0)})

where now

(∇g(x^{(0)})) = \frac{∂g_{i}}{∂x_{j}} |_{x=x^{(0)}}

As in the scalar case we base our iteration on the assumption that

l(x^{(k+1)}) = 0

Hence,

g(x^{(k)}) + ∇g(x^{(k)})(x^{(k+1)} − x^{(k)}) = 0

and given x^{(k)} we may determine the next iterate x^{(k+1)} by solving an n×n system of equations.

5.2.3. Unconstrained Minimization with Newton’s Method. In this chapter we are interested in minimizing functions of several variables. Analogously with the scalar variable case we may modify the above root finding method to determine maximum (or minima) of a function f(x_{1}, …, x_{n}). To compute an extreme point we require that ∇f = 0 hence we set

g(x) = (\frac{∂f}{∂x_{1}} , …, \frac{∂f}{∂x_{n}} )^{T}

Substituting

g_{i}(x) = \frac{∂f}{∂x_{i}}

into

g(x^{(k)}) + ∇g(x^{(k)})(x^{(k+1)} − x^{(k)}) = 0

produces

∇f(x^{(k)}) + Hf(x^{(k)})(x^{(k+1)} − x^{(k)}) = 0

where

(Hf(x^{(k)}))_{ij} = \frac{∂g_{i}}{∂x_{j}} = \frac{∂^{2}f}{∂x_{i}∂x_{j}}

Again we have a linear system for x^{(k+1)}.

5.3. Steepest Descent

Another form for Taylor’s formula in n-variables is given by

f(x + th) = f(x) + t∇f(x)h + \frac{t^{2}}{2} h^{T} Hf(x)h

where again (∇f(x))_{i} = ∂f/∂x_{i} and (Hf(x))_{ij} = ∂^{2}f/∂x_{i}∂x_{j}. Now t is a scalar and x + th is a ray emanating from the point x in the direction h. We can compute the derivative of the function f w.r.t. t as

\frac{df}{dt} (x + th) = ∇f(x)h + th^{T} Hf(x)f

Evaluating the derivative at the point t = 0 gives

\frac{df}{dt} (x) = ∇f(x)h

This quantity, known as the directional derivative, indicates how the function is changing at the point x in the direction h. Recall from calculus that the direction of maximum increase (decrease) of a function is in the direction of the gradient (negative gradient). This is readily seen from the formula for the directional derivative using the identity

∇f(x)h = ∥∇f(x)∥∥h∥ cos(θ)

We can assume without loss of generality that h is of unit length, i.e., ∥h∥ = 1. So the quantity on the right is a maximum when the vectors h and ∇f(x) point in the same direction so θ = 0.

This observation may be used to develop an algorithm for unconstrained function minimization. With an appropriate choice of the scalar step-size α, the iterations

x^{(k+1)} = x^{(k)} − α∇f(x^{(k)})

will converge (possibly slowly) to a minimum of the function f(x).