Recently, perusing the internet, I found an article which implements the trapezoidal rule in F#. 

open System

let main() =

//Function to integrate
let f x =
10.0*x*x
let trapezoidal a b N =
let mutable xi = a
let h = (b - a)/N
let mutable suma = h/2.0*(f(a)+f(b))
for x in 1 .. System.Convert.ToInt32(N) do
let mutable xi1 = xi + h
suma <- suma + h*f(xi1)
xi <- xi1
suma
//some usage example
let fromA = 0.0
let toB = System.Math.PI
let counter = 100.0
let result = trapezoidal fromA toB counter
Console.Write("Result = ")
Console.Write(result)
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. 

let h = (b - a)/N
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 .. System.Convert.ToInt32(N)|]
|> 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, 

float(i)*h+a
Then, 
suma + h*f(xi1)
will transform to 
acc + h*f(float(i)*h+a)
Where i is a current index,

So, the entire function will look as follows 

let trapezoidal2 a b N =
let h = (b - a)/N
let initialAccValue = h/2.0*(f(a)+f(b))
[|1 .. System.Convert.ToInt32(N)|]
|> 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.