Friday, September 4, 2015

Newton's Approximation

Way back in 1967 my school started an innovative new program. We could write programs for a computer at Imperial College in London, mail them in, have them run, with the results being returned to us the following week.

The biggest snag was that we didn't have a card-punch machine. So, instead we got cards that were semi-perforated--every second column (starting with 1) was perforated so that little rectangular holes could easily be punched out. Obviously, we had to know the EBCDIC code for all of the letters, numbers, etc.

The programs were written in Fortran so embedded spaces were meaningless. But if a line went on too long (a line was effectively limited to 37 characters) we needed a punch in column 6. We had special continuation cards since we couldn't create those ourselves.

Some of us decided to program a Newton's approximation of some easy function. A couple of the other kids chose cos(x) - x = 0. Because I was a smart-ass I chose cos(x) - sech(x) = 0.

I wanted to show an example of this for the class I'm teaching so, as usual, the first thing I did was to google "scala newton approximation" or something like that. One of the entries I got was this one:
But when I looked at that it was so un-Scala-like that I decided I would have to write my own. This is what I came up with:

/**
 * @author scalaprof
 * (c) Robin Hillyard (2015)
 */
case class Newton(f: Newtonian, guess: Double, maxTries: Int) {
    def solve: Either[String,Double] = solve(State(guess,maxTries))
    private def step(s: State): Either[State,Either[String,Double]] =
        f(s) match {
     case None => Right(Right(s.x))
     case Some(x) => s(x) match {
  case Left(e) => Right(Left(e))
  case Right(s1) => Left(s1)
     }
        }
import scala.annotation.tailrec
@tailrec private def solve(s: State): Either[String,Double] = step(s) match { case Right(r) => r case Left(l) => solve(l) } } case class State(x: Double, tries: Int) extends Function1[Double,Either[String,State]] { def apply(x: Double) = tries match { case 0 => Left("Failed to converge") case _ => Right(State(x,tries-1)) } } case class Newtonian(name: String, f: Double=>Double, dfbdx: Double=>Double, threshold: Double) extends Function1[State,Option[Double]] { def apply(s: State) = { val x = s.x val y = f(x) if (math.abs(y) > threshold) Some(x - y/dfbdx(x)) else None } } object Newton { def apply(f: Newtonian, start: Double): Newton = apply(f,start,100) def main(args: Array[String]): Unit = { val f = Newtonian("cos(x)-x (~1E-7)", {x => math.cos(x) - x},{x => -math.sin(x) - 1},1E-7) Newton(f,1.0).solve match { case Right(x) => println(s"""the solution to "${f.name}" is $x""") case Left(m) => println(s"error: $m") } } }
Notice that there are no variables in this program! All of the problem domain (including the definition of the function itself, its derivative and an appropriate tolerance value) is encapsulated in the Newtonian class and all of the logic of Newton's method is encapsulated in two classes: Newton which is the main driver and State which holds the current state of the run.

This makes good use of the Scala types Option and Either to cope with the logic aspects of the problem. There are really two places where this sort of logic is required: when running the actual approximation logic itself in method apply of class Newtonian: have we reached a point where we are close enough? Or should we return a new delta(x) value. Thus the return type is Option[Double] which can be either None or Some(delta).

The second place is in the apply method of the State class: if we've exhausted the maximum number of tries, we should return an error message, otherwise we should return a new value of State. This can be handled with an Either[String,State]. The values will thus be Left("Failed to converge") or Right(State(x,tries-1)).

In the solve method of class Newton, we sort all this logic out, returning an Either[String,Double] result.

Finally, we have the singular object Newton which includes the main program. This sets up the Newtonian for our particular problem, that's to say cos(x)-x=0 and tries to solve it. According to the result type, a message is printed on the console.

One other point to mention is that we do use recursion in this code. Programmers have been taught to avoid recursion because it is "less efficient than iteration." However, functional languages like Scala have a way of making recursion exactly as efficient as iteration, provided that the logic is expressed as a "tail" recursion. There is an annotation on the step method in Newton to ensure that our code really is tail-recursive.