Simpson Rule Sum in Functional Javascript

Photo by Stefan Grage on Unsplash

Both Riemann Sum and Trapezoid Sum has shown their respective capability of calculating the total area under certain curve. Their accuracy are only limited on how much you’re willing to put processing power into the calculation, the thinner the slices the more accurate the result should be. Another contender to challenge both of them is Simpson’s Rule, which in this article I’ll just write it as Simpson Sum.

Simpson Sum is still the method of exhaustion as implemented by other method, but there’s something profoundly different about this one. Instead of using linear geometric shape like rectangle or trapezoid, it uses rectangles with curved top. How? Let’s imagine that we have three different coordinate on a map like:

[
  {x: 3, y: 12},
  {x: 1, y:  5},
  {x: 5, y:  9}
]

Let’s assume that each of those coordinates follow a quadratic function where’s y = ax² + bx + c so each of that coordinate can be written as:

[
  {x: 3, y: 12}, // 12 = a*3² + b*3 + c
  {x: 1, y:  5}, //  5 = a*1² + b*1 + c
  {x: 5, y:  9}  //  9 = a*5² + b*5 + c
]

With linear programming, you could find the value of each a and b coefficient as well as the c constant. But as this article was intended to be short, just let me give you the result as {a: -1.25, b: 8.5, c: -2.25}. And if you draw it on Desmos, you’d find a graph like this.

Let me remind you that the equation we recently got may not be the only function that can touch each of those coordinates, but at least according to linear programming, that’s the closes approximation which can be derived. But let’s assume again that this is the equation we’re looking for, we can integrate the equation, set the limit from start to end, and get the appropriate result of integration. Since we don’t want to do it manually and looking for a way to let the computer does it, we require certain algorithm to help us finding the total area under a curve, specifically from a quadratic equation. Lucky for us, there’s a function that look like this:

It’s called the Simpson Integration Rule. The right hand of the function looks peculiar, how could sum-ing f(a) + 4f(m) + f(b) multiply it by the range (b - a), and divide them by 6 could be equal (or closely equal) to the real integration result? There’s a rigorous explanation here and here, but even with all of those help I can’t keep up well, so here I go trying to understand it intuitively:

square

Let’s say that we have a square that we can measure the area under y = 1 is square of 1. This does make sense, since both height and width are equally 1, thus the height * width must be square of 1. Now let’s try applying the Simpson Rule:

= ((2–1) / 6) * (f(1) + 4 * f(1.5) + f(2))
= ((2–1) / 6) * (  1  + 4 *    1   +  1  )
=   (1 / 6)   * (           6            )
= 1 // square of 1

Hmm,, it must be a coincidence right? Why would an algorithm which said to be specifically fit for a quadratic equation reach the right answer for area under horizontal line? Let’s have another try with a trapezoid like the graph below:

trapezoid

= ((2–1) / 6) * (f(1) + 4 * f(1.5) + f(2))
= ((2–1) / 6) * (  1  + 4 *   1.5  +  2  )
=   (1 / 6)   * (           9            )
= 1.5 // square of 1.5

The Simpson Rule got it right again!? I think it’s enough to prove that Simpson Rule does work in many circumstances other than quadratic equation alone. I may not understand the underlying math of how that function derived, but I’m assured that this function does work. But it can be seen that the former equation only use exactly 3 slice for any equation in any length. Could we do Simpson Rule with exhaustion method that use many more slices than just 3? Lucky again, there’s an equation just for that:

let isOdd = n => n % 2,
withAs = (obj, cb) => cb(obj),
makeArray = n => [...Array(n).keys()],
sum = arr => arr.reduce((r, i) => r + i),

simpsonSum = obj => withAs(
  (obj.end - obj.start) / obj.slice,
  width => width * sum(
    makeArray(obj.slice + 1)
    .map(i => [0, obj.slice]
      .includes(i) ? 1 : 2 * (1 + isOdd(i))
    ).map((i, j) => i * obj.func(
      obj.start + (width * j)
    ))
  ) / 3
)

simpsonSum({
  start: 0, end: 8, slice: 4,
  func: x => Math.pow(x, 1/2)
}) // get 14.8554

The function above can be called with more slices than just 3. As can be seen, the function call sample yields 14.8554 while the trapezoidSum from the last article gives us 14.5558 . It only proves that the function works just right. Now lets deal with a bunch of problems in the same manner of the last article:

problems = [
  { // a
    start: -1, end: 2, slice: 10,
    func: x => 3 * Math.pow(x, 2)
  }, // gets 9, while IC gets 9
  { // b
    start: 1, end: 8, slice: 10,
    func: x => Math.pow(x, 1/3)
  }, // gets 11.2496, while IC gets 11.25
  { // d
    start: 0, end: 2, slice: 10,
    func: x => 3*Math.pow(x, 3) - 2*x + 5
  }, // gets 18, while IC gets 18

  { // f
    start: Math.PI/6, end: Math.PI/2, slice: 10,
    func: x => Math.pow(Math.sin(x), 2) * Math.cos(x)
  }, // gets 0.2916, while IC gets 7/24 or 0.2916

  { // h
    start: 0, end: 1, slice: 10,
    func: x => x / Math.pow(Math.pow(x, 2) + 1, 2)
  }, // gets 0.2500, while IC gets 1/4 or 0.25

  { // i
    start: 0, end: 1, slice: 10,
    func: x => withAs(
      Math.pow(Math.E, x),
      exp => exp / (exp + 1)
    )
  }, // gets 0.6201, while IC gets ln(e+1)-ln(2) or 0.6201

  { // q
    start: 1, end: Math.E,
    slice: 10, func: Math.log
  }, // gets 0.9999, while IC gets 1

  { // r
    start: 1, end: Math.E, slice: 10,
    func: x => Math.pow(x, 3) * Math.log(x)
  }, // gets 10.2996, while IC gets 10.2996

  { // z
    start: 0, end: Math.PI/4, slice: 10,
    func: x => (
      1 + Math.pow(Math.sin(x), 2)
    ) / Math.pow(Math.cos(x), 2)
  } // gets 1.2146, while IC gets 2-pi/4 or 1.2146
]

problems.map(simpsonSum)

How impressive, with merely 10 slices each simpsonSum can achieve good accuracy as the last trapezoidSum . I think I’ll end the series here, looking forward on another math challenges to be solved with JS. Bye.