It is sometimes preferable to solve an integral numerically rather than analytically. An analytical solution is not always possible. Simpson's (Thomas Simpson, 1710 - 1761, English mathematician) rule is a method for approximating definite integrals. With it, an integral is approximated by a quadratic through three points within a given interval. Over the interval a...b, the integral of a function f(x) can be approximated by:
in which,∫ f(x)dx ≈ (y1 + 4·y2 + y3)·h/3
where,y1 = f(a) y2 = f(a + h) y3 = f(a + 2·h)
Suppose that we want to integrate over two equal intervals, a...m, m...b. Each interval would be divided into two segments of equal width h = (b - a)/4. Then, applying Simpson's rule to both intervals gives:h = (b - a)/2
∫ f(x)dx ≈ (y1 + 4·y2 + y3)·h/3 + (y3 + 4·y4 + y5)·h/3
≈ (y1 + 4·y2 + 2·y3 + 4·y4 + y5)·h/3
This idea can be extended to give a general version of the integration rule:
where,∫ f(x)dx ≈ (y1 + 4·y2 + 2·y3 + 4·y4 + ... + 2·y2·n-1 + 4·y2·n + y2·n+1)·h/3
While there are many known techniques for numerical integration, this one is easy to program, quickly executed, and very accurate for most practical cases. One implementation is shown below.h = (b - a)/(2·n) n = number of intervals in a...b
Zeno 1.2 is an interpreter for the Zeno programming language. It is an easy to learn and is suitable for educational purposes.
const PI : real := 2*arccos(0)
program
put "result = ", integrate( 0, PI/4, 1250 ):14:10
end program
% user defined function, returns value of integrand at x
function F( x : real ) : real
return cos(2*x)^(3/2)*sin(x)
end function
% compute the integral of F(x)
% a - lower limit, b - lower limit
% intervals - of simpson's rule integration
function integrate( a, b: real, intervals: int ): real
var sum, h, x: real
var n, i: int
assert intervals > 0 % valid argument?
n := 2*intervals % initialize
h := (b - a) / n
x := a
sum := F(a) + F(b) % first, last elements
for i := 2...n do % inner elements of summation
x := x + h % next x
if i mod 2 = 0 then % even index
sum := sum + 4*F(x)
else % odd index
sum := sum + 2*F(x)
end if
end for
return sum * h / 3 % multiply by common factor
end function
result = 0.1087094651