Refactoring F# Imperative Code Towards Declarative

Recently, perusing the internet, I found an article which implements the  trapezoidal rule in F#.
  1. open System  
  2.   
  3. let main() =  
  4.   
  5. //Function to integrate  
  6. let f x =  
  7. 10.0*x*x  
  8. let trapezoidal a b N =  
  9. let mutable xi = a  
  10. let h = (b - a)/N  
  11. let mutable suma = h/2.0*(f(a)+f(b))  
  12. for x in 1 .. System.Convert.ToInt32(N) do  
  13. let mutable xi1 = xi + h  
  14. suma <- suma + h*f(xi1)  
  15. xi <- xi1  
  16. suma  
  17. //some usage example  
  18. let fromA = 0.0  
  19. let toB = System.Math.PI  
  20. let counter = 100.0  
  21. let result = trapezoidal fromA toB counter  
  22. Console.Write("Result = ")  
  23. Console.Write(result)  
  24. main()  
What the code does is, over the course of 100 iterations it calculates a point, xi, and increments accumulator sums with the value of the function in this point multiplied by h.

Although there is nothing inherently wrong with this piece of code, its imperative fashion doesn't provide any benefit of functional approaches, such as expressiveness, minimizing mutability, reducing state side effects etc.

So, let's try to come up with trapezoidal2 which will take advantage of declarative programming.

First, let's extract initial values.
  1. let h = (b - a)/N  
  2. let initialAccValue = h/2.0*(f(a)+f(b))  
Notice how mutable variable sum was changed to non-mutable function. Also, as it's non-mutable, now, we've changed the name to the one that suggests that it's only the initial value for loop that can be replaced with the array, the contents of which, we'll iterate over using F# 's powerful collection traversing capabilities. As the goal is to return accumulator value the natural candidate is Array.fold,
  1. [|1 .. System.Convert.ToInt32(N)|]  
  2. |> Array.fold (fun acc i -> acc + //calculate increment here) initialAccValue  
The tricky part here is that at each step, we calculate not only accumulator value but also the value of a next point.
point(n) = point(n-1) + h

Luckily enough as h is constant we can also calculate it given first point and h

point(n) = point(0) + h*n

Or if we substitute our variables,
  1. float(i)*h+a  
Then,
  1. suma + h*f(xi1)  
will transform to
  1. acc + h*f(float(i)*h+a)  
Where i is a current index,

So, the entire function will look as follows
  1. let trapezoidal2 a b N =  
  2. let h = (b - a)/N  
  3. let initialAccValue = h/2.0*(f(a)+f(b))  
  4. [|1 .. System.Convert.ToInt32(N)|]  
  5. |> Array.fold (fun acc i -> acc + h*f(float(i)*h+a)) initialAccValue  
This is far more expressive and can be easily translated to other functional languages.