Class Notes: July 28, 2003 through August 1, 2003
CS 314

Optimization (continued from Friday's lecture)

• Assume we want to minimize f(x) and it has a unique minimum at p. where a <= p <= b and we start the search at p₀.

• If f'(p₀) > 0, p is to the left of p₀.

• If f'(p₀) < 0, p is to the right of p₀.

• If the derivative of f(x) is available, then we can use any of the methods we've covered to solve non-linear equations to solve f'(x) = 0.  (e.g. Regula-falsi, bisection, Newton Raphson, Secant, etc.)

Example:
f(x) = sin(x).  Find the minimum in the interval [-3pi/4, 0].  Use bisection on f'(x).
The answer is -1.5708

You could do numerical derivation if the derivative is hard to find.  Choose h carefully.  Any  method to solve f(x) = 0 can be used here.

Finding a Minimum for Multiple Variables

Steepest Descent or Gradient Method

Assume we want to minimize f(x) of n variables where x = (x₁, x₂, ..., x_n).  The gradient, ▼f(x), is a vector function defined as:
grad( f(x) ) = (df/dx₁, df/dx₂, ..., df/dx_n)
From the concept of gradient, we know that the gradient vector points in the directoin of the greatest increase of f(x).  Evaluate at a point to get the vector which points in the direction of greatest change.  Since we are looking for the minimum, we will use the negative of the gradient which will point in the direction of greatest decrease.

Basic Gradient Method
Start at point p₀ and move along the line in the direction ▼f(x) a small increment.  That is going to take you closer to the minimum.
    p₁ = p₀ - ▼f(p₀)*h        h is the small increment
    p_k = p_k-1 - ▼f(p_k-1)*h
When we reach p, the value of the gradient will be zero.

Example:
f(x) = x²+1
f'(x) = 2xi, where i is a unit vector in the x-direction

Let h = 0.1, x₀=-1
 

Numerical Solution of Differential Equations

Introduction

Example of a differential equation:
y' = k*y
dy/dt = k*y
dy = k*y dt
∫dy/y = ∫k dt
ln(y) = kt + k₂
y = e^{kt + k}₂
y = k₃e^kt

Very often, differential equations do not have an analytic solution so they have to be approximated using numerical methods.

Euler's Method

Let [a,b] be the interval which we want to find the solution of y' = f(y,t) with y(a) = y₀.  We will find a set of points: (t₀,y₀), ..., (t_n,y_n) that are used to approximate y(t) = y(t_k).
Also, t_k= a +k*h (where h is the step size) for k = {0,1,2,...,m} over [t₀, t₁, ..., t_m] with y(t₀) = 0.
Using Taylor's Expansion to approximate y(t) around t₀, we have:
        y(t) = y(t₀) + y'(t₀)*(t-t₀) + y''(c)*(t-t₀)²/2
We use this equation to obtain y(t₁):
        y(t₁) = y(t₀) + y'(t₀)*(t-t₀) + E(t)
If the step size is small enough, we can neglect the second order error.
        y₁ = y₀ + h*y'(t₀)
        y₁ = y₀ + h*f(t₀,y₀)
In general, t_k+1 = t_k +h
        y_k+1 = y_k + h*f(t_k,y_k) for k = 0, 1, 2, ..., m-1

Example:
y' = t²-y    y(0) = 1    h = 0.2

Heun's Method

We want to solve y'(t) = f( t, y(t) ) over [a,b] with y(t₀) = y₀.  We can integrate y'(t) over [t₀, t₁].
∫ y'(t) dt = ∫ f( t, y(t) ) dt
y(t₁) - y(t₀) = ∫ f( t, y(t) ) dt
Use numerical integration to approximate the integral in the right side.  Use trapezoidal rule.
y(t₁) - y(t₀) = (h/2)( f(t₀,y₀) + f(t₁,y₁) )
y(t₁) = y(t₀) + (h/2)( f(t₀,y₀) + f(t₁,y₁) )
Observe that we still need to know f(t₁,y(t₁)).  This involves y(t₁) that is what we want to solve.  To eliminate this circular reference, we use Euler's Approximation to approximate y(t₁)
        y(t₁) = y(t₀) + h*f(t₀,y₀)
The new equation is:
        y(t1) = y(t0) + (h/2)( f(t0,y0) + f(t1,y(t₀)+h*f(t₀,y₀) )
using y₀ = y(t₀) and y₁ = y(t₁)
    y₁ = y₀ + (h/2)( f(t₀,y₀) + f(t₀,y₀+h*f(t₀,y₀) )
    p_k+1 = y_k + h*f(t_k,y_k)
    y_k+1 = y_k + (h/2)( f(t_k,y_k) + f(t_k+1,p_k+1) )
Euler's Approximation is used as a predictor and the integral approximation is used as a correction.

Example: (same equation and initial values as the Euler's example)

Taylor Series Method

Using the Taylor Series to approximate the solution y(t₀+h) = y(t₀) + h*y'(t₀) + (h²/2!)(y''(t₀)) + ... + (hⁿ/n!)(y⁽ⁿ⁾(t₀))
We want to solve y' = f(t,y) when y(t₀) = a.  From the Taylor Expansion we can find the successive points y_k+1 = y_k + h*y'_k + (h²/2!)(y''_k) + ... + (hⁿ/n!)(y⁽ⁿ⁾_k).  We still need to compute y_k'' and y_k''', etc using y' = f(t,y).  The error in the approximation will be O(hⁿ⁺¹).  The higher the order of the Taylor Series, the smaller the error.

Example:
Solve y' = t²-y when y(0) = 1 and h = 0.2.  Use n=3.
y_k+1 = y_k + h*y'_k + (h²/2!)(y''_k) + (h³/3!)(y⁽³⁾_k)
y' = t²-y
y'' = 2t-y'
y''' = 2 - y''

y₁ = 1 + 0.2*-1 + 0.2²*1/2 + 0.2³*1/6 = 0.821333
y₁' = 0.2²-0.821333 = -0.781333
y₁'' = 2*0.2 - (-0.781333) = 1.181333
y₁''' = 2 - 1.181333 = 0.818667

y₂ = 0.821333 + 0.2(-0.781333) + 0.2²*1.181333/2 + 0.2³*0.818667/6 = 0.689785
y₂' = 0.4² - 0.689785 = -0.529785
y₂'' = 2*0.4 - (-0.529785) = 1.329785
y₂''' = 2 - 1.329785 = 0.670215

y₃ = 0.689785 + 0.2(-0.529785) + 0.2²*1.329875/2 + 0.2³*0.548683/6 = 0.611317
y₃' = 0.6² - 0.611317 = -0.251317
y₃'' = 2*0.6 - (-0.251317) = 1.451317
y₃''' = 2 - 1.451317 = 0.548683

y₄ = 0.611317 + 0.2(-0.251317) + 0.2²*1.451317/2 + 0.2³*0.548683/6 = 0.590812

Exact Solution:  y(0.8) = 0.5907
Taylor:  0.590812
Heun:  0.6001
Euler:  0.5123

To reduce approximation error, increase the order n of the expansion or reduce the value of h.  It is more effective to increase n because error is O(hⁿ⁺¹).

Runge-Kutta Method of Order 4

The Taylor Series Method gives a good approximation, but it requires the derivatives of the function.  The Runge-Kutta Method of Order 4 simulates the accuracy of the Taylor Series Method using an order of n=4 but it does not require the computation of derivatives.
y_k+1 = y_k + h*(f₁ + 2f₂ + 2f₃ + f₄)/6  where...
f₁ = f(t_k, y_k)
f₂ = f(t_k+h/2, y_k+f₁*h/2)
f₃ = f(t_k+h/2, y_k+f₂*h/2)
f₄ = f(t_k+h, y_k+f₃*h)

Example (same equation and initial conditions as previous example):
k = 0
y₀ = 1
f₁ = 0² - 1 = -1
f₂ = f(0.1, 0.9) = 0.1²-0.9 = -0.89
f₃ = f(0.1, 0.911) = 0.1² - 0.911 = -0.901
f₄ = f(0.2, 0.8198) = 0.2² - 0.8198 = -0.7798

k = 1
y₁ = 1 + 0.2*(-1 + 2*-0.89 + 2*-0.901 + -0.7798)/6 = 0.821273
f₁ = 0.2² - 0.821273 = -0.781273
f₂ = f(0.3, 0.743196) = 0.3² - 0.743196 = -0.653416
f₃ = f(0.3, 0.755958) = 0.3² - 0.755958 = -0.665958
f₄ = f(0.4, 0.688081) = 0.4² - 0.688081 = -0.528081

k = 2
y₂ = 0.821273 + 0.2*(-0.781273 + 2*-0.653146 + 2*-0.665958 + -0.528081)/6 = 0.689688

Systems of Differential Equations

Assume we have dx/dt = f(t,x,y), dy/dt = f(t,x,y), x(t₀) = x₀, and y(t₀) = y₀.  The solution to this system are functions x(t) and y(t) that when derivated and substituted in the system of equations give equality.

Example:
x' = x + 2y  x(0) = 6
y' = 3x + 2y  y(0) = 4
Solution:
x(t) = 4e^4t + 2e^-t
y(t) = 6e^4t - 2e^-t

Numerical Solution with Euler's Method

We can extend Euler's Method of a single differential equation to a system of equations.
dx/dt = f(t,x,y)  -->  dx = f(t,x,y) dt
dy/dt = f(t,x,y)  -->  dy = f(t,x,y) dt
dx_k = x_k+1 - x_k
dy_k = y_k+1 - y_k
dt_k = t_k+1 - t_k = h

x_k+1 - x_k = f(t_k, x_k, y_k)(t_k+1 - t_k)
y_k+1 - y_k = f(t_k, x_k, y_k)(t_k+1 - t_k)
x_k+1 = x_k + f(t_k, x_k, y_k)*h    x(t₀) = x₀
y_k+1 = y_k + f(t_k, x_k, y_k)*h    y(t₀) = y₀

The Euler Method does not give good accuracy because it approximates the Taylor Expansion only to the first derivative.

Runge-Kutta Method

The Runge-Kutta Method can be extended to systems of linear equations.  The formulas are as follows:

xk+1 = xk + (h/6)(f1 + 2f2 + 2f3 + f4)
yk+1 = yk + (h/6)(g1 + 2g2 + 2g3 + g4)

f1 = f(tk, xk, yk)
f2 = f(tk + h/2, xk + (h/2)f1, yk + (h/2)g1)
f3 = f(tk + h/2, xk + (h/2)f2, yk + (h/2)g2)
f4 = f(tk + h, xk + hf3, yk + hg3)

g1 = g(tk, xk, yk)
g2 = g(tk + h/2, xk + (h/2)f1, yk + (h/2)g1)
g3 = g(tk + h/2, xk + (h/2)f2, yk + (h/2)g2)
g4 = g(tk + h, xk + hf3, yk + hg3)

Higher Order Differential Equations

Higher order differential equations involve higher derivatives (e.g. x''(t) and y''(t)).  For example,

Example:
4x'' + 3x' + 5x = 2  x(0) = 1 and x'(0) = 3
x'' = (2 - 3x' - 5x)/4
y' = (2 - 3x' - 5x)/4  and  x' = y

Review for the Final Exam

20% first half
80% second half
Study: class notes, book, homework, do exercises from the book

Topics from the Second Half

• Least Squares Line
• Least Squares for non-linear equations
• Transformation for Data Linearization
• Polynomial Fitting
• Splines (the 5 properties, proof, how to use the formulas, end-point constraints)
• Numerical Differentiation (limit of difference quotient, central difference formula of both orders)
• Numerical Integration (Trapezoidal Rule and Simpson's Rule)
• Numerical Optimization (max & min, Golden Ratio Method, Minimization w/ Derivatives, Gradient Method)
• Solution of Differential Equations (Euler's Method, Heun's Method, Taylor Series, Runge-Kutta, Systems of DE, Higher Order DE)