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:

• A Scala implementation of Newton's Method for solving equations

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.

String Interpolation

Scala 2.10 introduced a nice feature for creating instances of types which can be defined using a String. It's called String Interpolation and its most common application is those strings that you see sprinkled through code like the following:

    throw new Exception(s"$x is too big for Int")

The expression beginning with "$" is replaced in the string by its value (in this case the evaluating x, whether that is a function or an identifier). Expressions including spaces require curly braces. If you follow the link referenced at the top, you will see how to customize string interpolation for your own purposes. If it is all completely clear to you by the time you get back here, then just ignore the rest of this blog post.

But if, like me, you find that comments like "A simple (buggy) implementation of this method could be:..." less than super-helpful, then stay right here and see a real example.

The example I'm going to present is the creation of a Rational number using the string form r"n/d" where n and d represent the numerator and denominator respectively. What Rational, I hear you say—there is no Rational object in Scala. Well, that's true (today at least). But we are going to suppose that we do have the following class (where most of the code in the middle has been replaced by ellipsis).

  case class Rational(n: Long, d: Long) extends Numeric[Rational] {
    ...  
  }

  object Rational {
    def apply(x: Long): Rational = new Rational(x,1)
    def apply(x: String): Rational = {
      val rRat = """^\s*(\d+)\s*(\/\s*(\d+)\s*)?$""".r
      x match {
        case rRat(n,_,d) => Rational(n.toLong,d.toLong)
        case rRat(n) => Rational(n.toLong)
        case _ => throw new Exception(s"invalid rational expression: $x")
      }
    }
    ...
  }

We will add the code described in the linked document to the object definition above. However, we will have to work the actual implementation a little better.

  object Rational {
    implicit class RationalHelper(val sc: StringContext) extends AnyVal {
      def r(args: Any*): Rational = {
        val strings = sc.parts.iterator
        val expressions = args.iterator
        val sb = new StringBuffer()
        while(strings.hasNext) {
          val s = strings.next
          if (s.isEmpty) {
            if(expressions.hasNext)
              sb.append(expressions.next)
            else
              throw new Exception("r: logic error: missing expression")
          }
          else
            sb.append(s)
        }
        if(expressions.hasNext)
          throw new Exception(s"r: logic error: ignored: ${expressions.next}")
        else
          Rational(sb.toString)
      }
   }
  ...
  }

You can probably eliminate the rather ugly logic errors because I don't think they can ever happen (although the mechanism which requires you to evaluate the next expression when you get an empty argument string is rather hokey and not well explained).

That's all there is to it. If you actually want to use it in another class file, then you will need to import Rational.RationalHelper.