27 December 2018

Numerical Integration - Simpson's Rule

Numerical integration is an operation that many older machines are not able to perform without add-on libraries (e.g. Math Pac for HP-41C) or without writing a program from scratch to do it. There are, of course, some exceptions to this; Hewlett Packard's HP-15C and HP-42S come straight to mind, as does the Casio fx-180P, for example, but these machines are actually the exception, not the rule.

There are many methods used to calculate definite integrals, that is integrals between known boundaries yielding a numerical result. One method that is a good trade-off between precision and complexity of application is Simpson's Rule, which attempts to find a quadratic of the form $a\cdot x^2+b\cdot x+c$ that fits or is close to our function. Finding the integral of this polynomial is very easy:
\[ \int_{x_1}^{x_2} \! (a\cdot x^2 + b\cdot x + c) \, \partial x = \frac a 3 (x_2^3-x_1^3) + \frac b 2 (x_2^2-x_1^2) + c\cdot (x_2 - x_1) \]
Basically, Simpson's Rule states:
\[ \int_a^{a+2h} \! f(x) \, \partial{x} \approx \frac{h}{3}(f(a)+4\cdot f(a+h)+f(a+2h)) \]
So, what does this mean? If we sample the function at the lower and upper boundaries and once more in the middle, we can get an approximation of the integral of the function by just adding the y-coordinates of the sample points and applying a few coefficients. Easy.

As with any approximation method that involves sampling, we can make this even more accurate by increasing the number of samples if we divide the region that we want to integrate into multiple smaller regions that we sample individually. This is thanks to (assuming a<c<b):
\[ \int_a^b \! f(x) \, \partial x = \int_a^c \! f(x) \, \partial x + \int_c^b \! f(x) \, \partial x \]
I wrote something for a few vintage machines that shows this method in action. The video below shows a Sharp PC-1360 and a Casio PB-410F running programs written in BASIC, followed by a SwissMicros DM41L and a Hewlett Packard HP-41CX running a similar program in RPN keystroke. You can get the listing for HP-/DM41 along with a .raw file and printable bar codes for use with a bar code wand here: simpson.zip

In this example, I'm getting all four machines to calculate:
\[ \int_0^1 \! (x^4-3x^3+x^2+1) \, \partial x \]
The exact answer to this is $ \frac{47}{60} \approx 0.7833333 $

In all examples I'm using 4 as the number of steps/divisions to use, which actually equates to $2^4=16$, i.e. I broke the [0;1] interval down into 16 smaller intervals (hence the progress counting up in increments of $\frac 1 {16} $= 6.25%).

No comments: