Created composite Simpson's integration method.Tests included

Author: ggkogkou

Maths/SimpsonIntegration.js

@@ -0,0 +1,69 @@
+/*
+*
+* @file
+* @title Composite Simpson's rule for definite integral evaluation
+* @author: [ggkogkou](https://github.com/ggkogkou)
+* @brief Calculate definite integrals using composite Simpson's numerical method
+*
+* @details The idea is to split the interval in an EVEN number N of intervals and use as interpolation points the xi
+* for which it applies that xi = x0 + i*h, where h is a step defined as h = (b-a)/N where a and b are the
+* first and last points of the interval of the integration [a, b].
+*
+* We create a table of the xi and their corresponding f(xi) values and we evaluate the integral by the formula:
+* I = h/3 * {f(x0) + 4*f(x1) + 2*f(x2) + ... + 2*f(xN-2) + 4*f(xN-1) + f(xN)}
+*
+* That means that the first and last indexed i f(xi) are multiplied by 1,
+* the odd indexed f(xi) by 4 and the even by 2.
+*
+* N must be even number and a<b. By increasing N, we also increase precision
+*
+* More info: [Wikipedia link](https://en.wikipedia.org/wiki/Simpson%27s_rule#Composite_Simpson's_rule)
+*
+*/
+
+
+function integralEvaluation(N, a, b, func) {
+    // Check if N is an even integer
+    let isNEven = (N) => {
+        if(N%2  ===  0)
+            return true;
+        return false;
+    }
+    if(!Number.isInteger(N) || Number.isNaN(a) || Number.isNaN(b)) { throw new TypeError("Expected integer N and finite a, b"); }
+    if(!isNEven) { throw "N is not an even number"; }
+
+    // Check if a < b
+    if(a > b) { throw "a must be less or equal than b"; }
+    if(a === b) return 0;
+
+    // Calculate the step h 
+    const h = (b - a) / N;
+
+    // Find interpolation points
+    let xi = a; // initialize xi = x0
+    let pointsArray = [];
+
+    // Find the sum {f(x0) + 4*f(x1) + 2*f(x2) + ... + 2*f(xN-2) + 4*f(xN-1) + f(xN)}
+    let temp;
+    for(let i = 0; i < N+1; i++){
+        if(i === 0 || i === N) temp = func(xi);
+        else if(i%2 === 0) temp = 2*func(xi);
+        else temp = 4*func(xi);
+
+        pointsArray.push(temp);
+        xi += h;
+    }
+
+    // Calculate the integral
+    let result = h/3;
+    temp = 0;
+    for(let i=0; i < pointsArray.length; i++) temp += pointsArray[i];
+
+    result *= temp;
+
+    if (Number.isNaN(result)) { throw "Result is NaN. The input interval doesn't belong to the function's domain"; }
+
+    return result;
+}
+
+export { integralEvaluation }

Maths/test/SimpsonIntegration.test.js

@@ -0,0 +1,16 @@
+import { integralEvaluation } from "./SimpsonIntegration"; 
+
+test("Should return the integral of f(x) = sqrt(x) in [1, 3] to be equal 2.797434", () => {
+    const result = integralEvaluation(16, 1, 3, (x) => {return Math.sqrt(x);});
+    expect(Number(result.toPrecision(7))).toBe(2.797434);
+});
+
+test("Should return the integral of f(x) = sqrt(x) + x^2 in [1, 3] to be equal 11.46410161", () => {
+    const result = integralEvaluation(64, 1, 3, (x) => {return Math.sqrt(x) + Math.pow(x, 2);});
+    expect(Number(result.toPrecision(10))).toBe(11.46410161);
+});
+
+test("Should return the integral of f(x) = log(x) + Pi*x^3 in [5, 12] to be equal 15809.9141543", () => {
+    const result = integralEvaluation(128, 5, 12, (x) => {return Math.log(x) + Math.PI*Math.pow(x, 3);});
+    expect(Number(result.toPrecision(12))).toBe(15809.9141543);
+});