Elemar Júnior

Posted on

# How to parse, simplify, differentiate and evaluate/solve (using Newton's method) equations and expressions in F#.

This was orginally posted in my blog

# The Object Model

The first step to parse expressions and equations is to create a "good-enough" object model.

type Expr =
| X
| Const of value: double
| Add   of Expr * Expr
| Sub   of Expr * Expr
| Mult  of Expr * Expr
| Div   of Expr * Expr
| Pow   of Expr * Expr
| Neg   of Expr

// 2 + 2 / 2
let sample = Add(Const(2.), Div(Const(2.), Const(2.))


It is fantastic how easy it is when using F#.

# Parsing

In this implementation, I use F# Active Patterns to do the parsing.

open System

let (|Digit|_|) = function
| x::xs when Char.IsDigit(x) ->
Some(Char.GetNumericValue(x), xs)
| _ -> None

let (|IntegerPart|_|) input =
match input with
| Digit(h, t) ->
let rec loop acc = function
| Digit(x, xs) -> loop ((acc * 10.) + x) xs
| xs -> Some(acc, xs)
loop 0. input
| _ -> None

"10" |> List.ofSeq |> (|IntegerPart|_|)

let (|FractionalPart|_|) = function
| '.'::t->
let rec loop acc d = function
| Digit(x, xs) -> loop ((acc * 10.) + x) (d * 10.) xs
| xs -> (acc/d, xs)
Some(loop 0. 1. t)
| _ -> None

"10" |> List.ofSeq |> (|FractionalPart|_|)
".34" |> List.ofSeq |> (|FractionalPart|_|)

let (|Number|_|) = function
| IntegerPart(i, FractionalPart(f, t)) -> Some(i+f, t)
| IntegerPart(i, t) -> Some(i, t)
| FractionalPart(f, t) -> Some(f, t)
| _ -> None

"10" |> List.ofSeq |> (|Number|_|)
".35" |> List.ofSeq |> (|Number|_|)
"10.35" |> List.ofSeq |> (|Number|_|)

let parse (expression) =
let rec (|Expre|_|) = function
| Multi(e, t) ->
let rec loop l = function
| '+'::Multi(r, t) -> loop (Add(l, r)) t
| '-'::Multi(r, t) -> loop (Sub(l, r)) t
| [] -> Some(l, [])
| _ -> None
loop e t
| _ -> None
and (|Multi|_|) = function
| Atom(l, '*'::Powi(r, t)) -> Some(Mult(l, r), t)
| Atom(l, '/'::Powi(r, t)) -> Some(Div(l, r), t)
| Powi(e, t) -> Some(e, t)
| _ -> None
and (|Powi|_|) = function
| '+'::Atom(e, t) -> Some(e, t)
| '-'::Atom(e, t) -> Some(Neg(e), t)
| Atom(l, '^'::Powi(r, t)) -> Some(Pow(l, r), t)
| Atom(e, t) -> Some(e, t)
| _ -> None
and (|Atom|_|) = function
| 'x'::t -> Some(X, t)
| Number(e, t) -> Some(Const(e), t)
| '('::Expre(e, ')'::t) -> Some(e, t)
| _ -> None

let parsed = (expression |> List.ofSeq |> (|Expre|_|))

match parsed with
| Some(result, _) -> result
| None -> failwith "failed to parse expression"

exec 0. (parse "2+2") // 4
exec 2. (parse "x^3")
parse "x^2-3" // Sub(Pow(X, Const(2.)), Const(3.)


# Simplifying and Evaluating

The following code can simplify equations/expressions removing steps to solve it.

let rec simplify e =
let result =
match e with
| Add(Const(0.), r) -> simplify r
| Add(l, Const(0.)) -> simplify l
| Add(Const(l), Const(r)) -> Const (l + r)
// sub
| Sub(Const(0.), r) -> Neg (simplify r)
| Sub(l, Const(0.)) -> l
| Sub(Const(l), Const(r)) -> Const (l - r)
| Sub(X, r) -> Sub (X, simplify r)
| Sub(l, X) -> Sub (simplify l, X)
| Sub(l, r) -> (Sub(simplify l, simplify r))
// mult
| Mult(Const(0.), _) -> Const(0.)
| Mult(_, Const(0.)) -> Const(0.)
| Mult(Const(1.), r) -> r
| Mult(l, Const(1.)) -> l
| Mult(Const(l), Const(r)) -> Const (l * r)
| Mult(l, r) when l = r -> (Pow (simplify l, Const(2.)))
| Mult(Pow(b, p), r) when b = r -> Pow (simplify b, (simplify (Add((simplify p), Const(1.)))))
| Mult(X, r) -> Mult (X, simplify r)
| Mult(l, X) -> Mult (simplify l, X)
| Mult(l, r) -> (Mult(simplify l, simplify r))
// div
| Div(Const(0.), _) -> Const(0.)
| Div(l, Const(1.)) -> l
| Div(Const(l), Const(r)) -> Const (l / r)
| Div(X, r) -> Div (X, simplify r)
| Div(l, X) -> Div (simplify l, X)
| Div(l, r) -> simplify (Div(simplify l, simplify r))
// pow
| Pow(_, Const(0.)) -> Const(1.)
| Pow(b, Const(1.)) -> simplify b
| Pow(Const(l), Const(r)) -> Const(System.Math.Pow(l, r))
| Pow(X, r) -> Pow (X, simplify r)
| Pow(l, X) -> Pow (simplify l, X)
| Pow(b, p) -> (Pow(simplify b, simplify p))
// neg
| Neg(Const(k)) -> Const (-k)
| Neg(X) -> Neg(X)
| Neg(x) -> (Neg(simplify x))
//
| other -> other

if (result = e)
then result
else simplify result

simplify (Mult(Mult(X, X), X))
simplify (Pow(Const(2.), Const(3.)))
simplify (Mult(Const(2.), X))

I love local functions! The simplification process works as an evaluator for expressions. With equations, the process will stop when there are no more possible simplification steps to take.

let exec x expr =
let rec replaceX = function
| Sub(l, r)  -> Sub(replaceX l, replaceX r)
| Mult(l, r) -> Mult(replaceX l, replaceX r)
| Div(l, r)  -> Div(replaceX l, replaceX r)
| Pow(l, r)  -> Pow(replaceX l, replaceX r)
| Neg(e)     -> Neg(replaceX e)
| Const(v)   -> Const(v)
| X          -> Const(x)

match simplify (replaceX expr) with
| Const(result) -> result
| _ -> failwith "impossible to execute"

// resulta 8
Pow(Const(2.), X) |> exec 3.


# Differentiating

Newton's method will need derivatives to work. So, let's produce it.

let rec deriv = function
| X          -> Const(1.)
| Const(_)   -> Const(0.)
| Sub(l, r)  -> Sub(deriv l, deriv r)
| Mult(l, r) -> Add(Mult(deriv l, r), Mult(l, deriv r))
| Neg(v)     -> Neg(deriv v)
| Pow(b, e)  -> Mult(e, simplify (Pow(b, Sub(e, Const(1.)))))
| _          -> failwith "expression not supported."

deriv (Pow(X, Const(3.)))


# Welcome, Mr. Newton

We now already have all the elements we need to solve equations using Newton's method.

Here is my implementation.

let newton expr guess error maxdepth =
let o = parse expr
let d = deriv o
let eq = Sub(X, Div(o, d))

let rec iter g depth =
if depth > maxdepth
then failwith "maxdepth exceeded."
else
let newg = exec g eq
printfn "%A" g

if (abs (newg - g) < error)
then newg
else iter newg (depth + 1)

iter guess 0

newton "x^3-27" 5. 0.000001 100 // 3


The parameters are the equation we need to solve, a solution guess, an acceptable error, and iterations.

We can make it simpler to use:

let solve expr =
newton expr 100. 0.00001 100

solve "x^2-9"        // 3
solve "3*x^2-4*x+1"  // 1


That's it

This post was written just for fun. I did it trying to learn F#. Do you like? Could I do something different? Please, give me your feedback.