The Science of Programming/Peaks and Valleys

According to CME, To find a minima or maxima of a curve with independent variable x, you take the derivative, set the derivative to zero, and then solve for x. At that point, the slope of the original is zero and therefore must be the highest point on the peak or the lowest point in a valley. As CME states, you don't know for sure if the zero slope represents a peak or if it represents a valley, but it must be one or the other.

The process for solving for x is much like the process for simplifying an expression; it is a complex task beyond the scope of this text. So we will not be able to find minima/maxima symbolically. However, thanks to Sir Isaac Newton, we will be able to perform this task numerically.

For the rest of this chapter, we will dispense with our object approach and take a functional approach, as the whole purpose for using objects was to perform symbolic tasks. Adapting our previous object approach to perform numeric tasks such as finding minima/maxima is left as an exercise.

A Blast from the Past

To begin, we must learn how to represent polynomials and their derivatives as numeric Sway functions. Remember our ratio function from Chapters 3 and 4? We used it to compute ${\frac {dy}{dx}}$ numerically:

   function ratio(x,dx,y)
       {
       var dy = y(x + dx) - y(x);
       dy / dx;
       }

In Sway, not only are we allowed to return numbers, strings and environments (objects) from functions, we can also return functions themselves. We will use this fact to build numeric derivative functions.

---

At this point, you should read about returning local functions in the Sway Reference Manual.

---

So rather than use the ratio function which requires us to pass in the same y function over and over, we will define a function to return a customized ratio that takes only the x value as an argument. Looking at the code will help immensely:

   function derivative(y,dx)
       {
       function ratio(x)
           {
           var dy = y(x + dx) - y(x);
           dy / dx;
           }
       }

Since the last expression in the derivative function is the definition of ratio, it is the ratio function that is the return value for derivative. Now, we can pass and x value to both the original function and its derivative.

As always, let's test, using this polynomial as a guinea pig:

   $3x^{2}-10x+5$

Representing this polynomial numerically as a Sway function, we have:

   function y(x)
       {
       3 * (x ^ 2) - (10 * x) + 5;
       }
     
   var y' = derivative(y,0.000001);

We would expect the symbolic derivative of y to be:

   $6x-10$

Therefore, we would expect y(2) to yield -3 and y'(2) to yield 2.

   sway> y(2);
   INTEGER: -3
   sway> y'(2);
   REAL_NUMBER: 2.0000029988

The imprecision in the second result is again due to working with real numbers rather than integers.

A good fixation

A necessary next step in finding minima/maxima numerically is to be able to find a fixed point of a function. If a function has a fixed point (and not all functions do), then there is a value x such that:

   $f(x)=x$

For example, the cosine function has a fixed point near 0.7390851332:

   sway> cos(0.7390851332);
   REAL_NUMBER: 0.7390851332

I say near because Sway does not show all the significant digits, by default.

Finding a fixed point is algorithmically simple:

   function fixedPoint(f,x)
       {
       if (f(x) == x)
           {
           x;
           }
       else
           {
           fixedPoint(f,x + f(x) / 2);
           }
       }

If, for the given x, $f(x)!=x$, we try again with a new value of x. Conveniently, the average of x and $f(x)$ is a mathematically sound value for our next try.

Although this code illustrates the approach well, there are a few things wrong with it. The first problem is that f(x) is calculated twice. For efficiency's sake, we should calculate it once:

   function fixedPoint(f,x)
       {
       var next = f(x);
       if (next == x)
           {
           x;
           }
       else
           {
           fixedPoint(f,x + next / 2);
           }
       }

The next problem is due to the limited precision of real numbers in Sway. Because real numbers cannot be represented exactly, it is usually unwise to make equality comparisons between real numbers at the extremes of their precision. Therefore, unless you know better (and in the case of fixedPoint, you don't know better), you should check instead if the numbers are within some small amount:

   var FIXED_POINT_DELTA = 1e-10;

   function fixedPoint(f,x)
       {
       var next = f(x);
       if (abs(x - next) < FIXED_POINT_DELTA)
           {
           x;
           }
       else
           {
           fixedPoint(f,x + next / 2);
           }
       }

The abs (absolute value) function is used since x - next might be a large negative value, which would also be less than $1*10^{-10}$.

We need an initial value for x when we call fixedPoint. It turns out that 1.0 is almost always a good guess, so let's hard-wire it. Here is a neat trick for doing that:

   function fixedPoint(f)
       {
       function helper(f,x)
           {
           var next = f(x);
           if (abs(x - next) < FIXED_POINT_DELTA)
               {
               x;
               }
           else
               {
               helper(f,x + next / 2);
               }
           }
       helper(f,1.0);
       }

The essence of the trick is to rename your function for which you wish to fix some arguments to a nice name like helper. You then wrap that function with a function having the original name, making helper a local function. Then the last thing you do in the wrapper function is call the helper, fixing the values of the desired argument.

We can improve this latest version by noting that fixedPoint's formal parameter f is bound to the same value as helper's formal parameter f, so we can remove helper's version in both the definition and call:

   function fixedPoint(f)
       {
       function helper(x)
           {
           var next = f(x);
           if (abs(x - next) < FIXED_POINT_DELTA)
               {
               x;
               }
           else
               {
               helper(x + next / 2);
               }
           }
       helper(1.0);
       }

Don't forget to change the recursive calls as well.

Newton's Transformer

Newton came up with a clever way to find out where a polynomials (and other differentiable functions) have a value of zero. He said that a zero of some function f is also a fixed point of a new function derived from f. This new function, $f_{NT}$, is:

   $f_{NT}(x)=x-{\frac {f(x)}{f\,'(x)}}$

Thus, to find a zero of some function, we generate a new function using Newton's Transform and pass the transformed function to our fixed point finder. The value it returns is our zero!

First, let's implement Newton's Transform:

   function NewtonsTransform(f)
       {
       var f' = derivative(f,0.000001);
       function y(x)
           {
           x - (f(x) / f'(x));
           }
       }

Now let's test our polynomial:

   $y=3x^{2}-10x+5$

We implement it using a Sway function, as before, and create the derivative function:

   function y(x)
       {
       3 * (x ^ 2) - (10 * x) + 5;
       }
     
   var y' = derivative(y,0.000001);

Where does function y produce a zero?

   sway> fixedPoint(NewtonsTransform(y));
   REAL_NUMBER: 0.6125741133
   sway> y(0.6125741133);
   REAL_NUMBER: -1.441567e-10

Given our inherent imprecision, 0.6125741133 does seem to be near a zero for y. How about y'?

   sway> fixedPoint(NewtonsTransform(y));
   REAL_NUMBER: 1.6666661669
   sway> y(1.6666661669);
   REAL_NUMBER: 1.776357e-09

Again, close enough to zero for government work.

Questions

1. Modify the symbolic system to handle numeric derivatives. To make numeric derivatives compatible, the derivative function should return an object with toString and value methods.