|
| 1 | +/* |
| 2 | +* |
| 3 | +* @file |
| 4 | +* @title Composite Simpson's rule for definite integral evaluation |
| 5 | +* @author: [ggkogkou](https://github.com/ggkogkou) |
| 6 | +* @brief Calculate definite integrals using composite Simpson's numerical method |
| 7 | +* |
| 8 | +* @details The idea is to split the interval in an EVEN number N of intervals and use as interpolation points the xi |
| 9 | +* 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 |
| 10 | +* first and last points of the interval of the integration [a, b]. |
| 11 | +* |
| 12 | +* We create a table of the xi and their corresponding f(xi) values and we evaluate the integral by the formula: |
| 13 | +* I = h/3 * {f(x0) + 4*f(x1) + 2*f(x2) + ... + 2*f(xN-2) + 4*f(xN-1) + f(xN)} |
| 14 | +* |
| 15 | +* That means that the first and last indexed i f(xi) are multiplied by 1, |
| 16 | +* the odd indexed f(xi) by 4 and the even by 2. |
| 17 | +* |
| 18 | +* N must be even number and a<b. By increasing N, we also increase precision |
| 19 | +* |
| 20 | +* More info: [Wikipedia link](https://en.wikipedia.org/wiki/Simpson%27s_rule#Composite_Simpson's_rule) |
| 21 | +* |
| 22 | +*/ |
| 23 | + |
| 24 | +function integralEvaluation (N, a, b, func) { |
| 25 | + // Check if N is an even integer |
| 26 | + let isNEven = true |
| 27 | + if (N % 2 !== 0) isNEven = false |
| 28 | + |
| 29 | + if (!Number.isInteger(N) || Number.isNaN(a) || Number.isNaN(b)) { throw new TypeError('Expected integer N and finite a, b') } |
| 30 | + if (!isNEven) { throw Error('N is not an even number') } |
| 31 | + if (N <= 0) { throw Error('N has to be >= 2') } |
| 32 | + |
| 33 | + // Check if a < b |
| 34 | + if (a > b) { throw Error('a must be less or equal than b') } |
| 35 | + if (a === b) return 0 |
| 36 | + |
| 37 | + // Calculate the step h |
| 38 | + const h = (b - a) / N |
| 39 | + |
| 40 | + // Find interpolation points |
| 41 | + let xi = a // initialize xi = x0 |
| 42 | + const pointsArray = [] |
| 43 | + |
| 44 | + // Find the sum {f(x0) + 4*f(x1) + 2*f(x2) + ... + 2*f(xN-2) + 4*f(xN-1) + f(xN)} |
| 45 | + let temp |
| 46 | + for (let i = 0; i < N + 1; i++) { |
| 47 | + if (i === 0 || i === N) temp = func(xi) |
| 48 | + else if (i % 2 === 0) temp = 2 * func(xi) |
| 49 | + else temp = 4 * func(xi) |
| 50 | + |
| 51 | + pointsArray.push(temp) |
| 52 | + xi += h |
| 53 | + } |
| 54 | + |
| 55 | + // Calculate the integral |
| 56 | + let result = h / 3 |
| 57 | + temp = 0 |
| 58 | + for (let i = 0; i < pointsArray.length; i++) temp += pointsArray[i] |
| 59 | + |
| 60 | + result *= temp |
| 61 | + |
| 62 | + if (Number.isNaN(result)) { throw Error('Result is NaN. The input interval doesnt belong to the functions domain') } |
| 63 | + |
| 64 | + return result |
| 65 | +} |
| 66 | + |
| 67 | +export { integralEvaluation } |
0 commit comments