Let’s Integrate – Using APL
Vector 8:1, page 22
Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 22
Introduction¶
The great beauty about APL, and hence its great value to teachers, is the ability to do such a lot with a few small programs. A few simple tools provided for eager students to use can give as much stimulus as some very exotic teaching software. The express aim of this article is to do just that, choosing the topic of integration. No matter how rusty our Mathematics is, most of us can remember what integration is all about - ‘areas under curves’ is an adequate description for most purposes.
Where to start? Well, since everything is about areas under curves, and since the simplest area to calculate is that of a rectangle, a useful tool is to be able to specify a set of contiguous intervals (each interval of which will support a rectangle of known height).
A GRID Function¶
As in one of my previous articles, design is quite important. I require here a function, call it GRID, which will calculate the end-points of N sub-intervals of a specified interval [a,b].
∇ R←GRID A_B_N;A;B;N;⎕IO
[1] ⎕IO←1 ⋄ A←A_B_N[1] ⋄ B←A_B_N[2] ⋄ N←A_B_N[3]
[2] R←A+(÷N)×(B-A)×0,⍳N
∇
GRID 0 10 10
0 1 2 3 4 5 6 7 8 9 10
GRID 0 2 20
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6
1.7 1.8 1.9 2
Note that the last of the three arguments supplies the number of intervals - not the number of points. Of course this simple function could be made more robust if required - obvious design features might include:
(a) testing that the right argument is numeric
(b) testing that the right argument has three elements
(c) testing that A≤B and that N is a positive integer.
Although this aspect is not central to the theme of this article, some readers may find it instructive to see just what is involved in providing such error traps. Perhaps the easiest way to do this is to use ‘execute’ as I suggested in my previous article. This is how GRID might be expanded to include the above suggestions.
∇ R←GRID A_B_N;A;B;N;⎕IO
[1] A_B_N←,A_B_N ⍝ force right argument to be a vector
[2] ⍎(0≠1↑0⍴A_B_N)/'''Right Argument must be numeric''⋄→0'
[3] ⍎(3≠⍴A_B_N)/'''Right Argument must have 3 elements''⋄→0'
[4] ⎕IO←1 ⋄ A←A_B_N[1] ⋄ B←A_B_N[2] ⋄ N←A_B_N[3]
[5] ⍎(A>B)/'''Interval mis-specified''⋄→0'
[6] ⍎((N≠⌊N)∨N≤0)/'''Number of Intervals must be a positive integer''⋄→0'
[7] R←A+(÷N)×(B-A)×0,⍳N
∇
Having re-written the function with error trapping, it all has to be tested!
GRID 'FRED'
Right Argument must be numeric
GRID 2 2⍴1 2 3 4
Right Argument must have 3 elements
GRID 2 1 4
Interval mis-specified
GRID 1 2 0
Number of Intervals must be a positive integer
GRID 1 2 10.5
Number of Intervals must be a positive integer
GRID 0 2 4
0 0.5 1 1.5 2
(Readers who have not gone through this exercise may now realise why it is so easy to crash software. There is no particular difficulty in error-trapping individual functions, but some tests could be easily left out in what is clearly a laborious process.)
Applying GRID¶
Back to our main theme - integration. As an illustration of how effective just this one tool is, let’s take a simple function to integrate - say the function FUN written below, given by the expression x exp(−x), or in APL, X×*-X, and integrate it from 0 to 2
As it happens, my integration is just about up to this calculation (making numerical integration redundant, which however is not the point) and the answer is given by [−x exp(−x) −exp(−x)] evaluated at 2 and 0. Hence the correct answer is 1 − 3 exp(−2)
As a simple approximation to the required integral, the range of integration may be divided into ten intervals and the function evaluated at all these points:
GRID 0 2 10
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
FUN GRID 0 2 10
0 0.1637461506 0.2681280184 0.3292869817 0.3594631713
0.3678794412 0.3614330543 0.3452357495 0.3230344288
0.2975379988 0.2706705665
giving the heights of the rectangles. Hence an approximation is the sum of the heights times the interval widths:
Trying finer grids will of course give better approximations:-
Given knowledge about the behaviour of the function (increasing from 0 to 1 and then decreasing from 1 onwards) grids can be chosen (always to include the turning point x=1 as an end-point) so that upper and lower bounds can be obtained simply by checking which end-point of each interval has the smaller (larger) function value. The easiest way to do this is to turn the result of GRID into a two column matrix, each row specifying the end-points of each interval.
⎕PP←4
MAT←GRID 0 2 10
MAT←(¯1↓MAT),[1.5]1↓MAT
MAT one row for each interval
0 0.2
0.2 0.4
0.4 0.6
0.6 0.8
0.8 1
1 1.2
1.2 1.4
1.4 1.6
1.6 1.8
1.8 2
FUN MAT
0 0.1637
0.1637 0.2681
0.2681 0.3293
0.3293 0.3595
0.3595 0.3679
0.3679 0.3614
0.3614 0.3452
0.3452 0.323
0.323 0.2975
0.2975 0.2707
.2×+/⌈/FUN MAT Calculate Maximum of each row and use
0.6367
.2×+/⌊/FUN MAT Calculate Minimum of each row and use
0.5437
Hence the integral lies between .5437 and .6367. Now to close in on it by repeating the operation with a finer grid. But first, note that really I have been rather naughty in applying FUN after creating MAT rather than creating MAT after using FUN on the result of GRID. So, subsequently, I shall reverse the order of operations. (To do so becomes important if the function to be integrated is complex and evaluations take a long time.)
Hence the integral lies between .5893 and .5986. Once more, take a finer grid:
So with one useful function and a few APL array operations, it is possible to explore the very foundations of integration, whereby the existence of the integral is proven by a pincer movement from above and below. (There are some bizarre mathematical functions for which this process fails!)
The Trapezium Rule¶
After the theory of integration, then consider some more simple techniques of numeric integration. The first method used above, merely approximated the function using a step-function - i.e. lots of horizontal bars. Since most functions integrated in practice are smooth, it seems reasonable to approximate by connecting the function points on the grid with straight lines. This leads to using the trapezium rule to evaluate the area of each trapezium standing on each interval. Each such trapezium has area equal to the average of the heights times the width. If this is applied to all of them, then a formula for the integral is:
With a view to the next method, rather than taking the sum of all the f values and subtracting .5 times the sum of the first and last values, we will generate a vector of coefficients .5 1 1 1 1 ,...,.5 using the expression .5,((N-1)⍴1),.5. This can then be multiplied and added to the f values using +.×:
This is the basis of the trapezium rule, which can now be incorporated into a function. To do so raises the interesting question of how to specify the function to be integrated. The only way possible with first-generation APL interpreters is to pass in the name of the function as a string variable. Consequently, the execute function (see Education Vector April 1991) has to be used.
There is however a problem with passing function names into an APL function, namely that care must be taken to make sure that no local variable names can coincide with the function name passed in. (If such a clash does occur, then what is executed is the LOCAL variable.) A common way of avoiding such a pitfall is to use variable names using ∆ or ⍙. So with this in mind, our integration function will have as its left argument the name of the function, and as its right argument the parameter of the GRID function to be used.
Here is the function TR_INT. The key part is line 3 to the right of ⍎, where a string vector is constructed consisting of the function name, followed by a space (very important) and the GRID command, all within quotes.
∇ R←∆FN TR_INT ∆A_B_N;∆A;∆B;∆N;∆WT
[1] ∆A←∆A_B_N[1] ⋄ ∆B←∆A_B_N[2] ⋄ ∆N←∆A_B_N[3]
[2] ∆WT←0.5,((∆N-1)⍴1),0.5
[3] R←(÷∆N)×(∆B-∆A)×∆WT+.×⍎∆FN,' GRID ∆A_B_N'
∇
'FUN' TR_INT 0 2 100
0.593956
'FUN' TR_INT 0 2 1000
0.593994
Simpson’s Rule¶
So far so good: one top level function using a suitable tool. What next? Well, as the trapezium rule supersedes the simple step-function approximation, so Simpson’s rule which approximates a curve by fitting quadratic curves to each interval normally supersedes the trapezium rule. The coefficients multiplying the end-points of an assumed even number of intervals are now 1 4 2 4 2 4 2….4 1, divided by 3 - otherwise the function is very similar - I have called it SMP_INT.
∇ R←∆FN SMP_INT ∆A_B_N;∆A;∆B;∆N;∆WT
[1] ∆A←∆A_B_N[1] ⋄ ∆B←∆A_B_N[2] ⋄ ∆N←∆A_B_N[3]
[2] ∆WT←1,((∆N-1)⍴ 4 2),1
[3] R←(÷3×∆N)×(∆B-∆A)×∆WT+.×⍎∆FN,' GRID ∆A_B_N'
∇
'FUN' SMP_INT 0 2 100
0.593994
Note that this is accurate to the number of decimal places shown! Note also that in some cases it is possible to evaluate the integral required by supplying a definition, for example:
Here is another example using the circular functions. First of all, multiples of pi are given by ○ right argument. Then 1○ right argument specifies the sine function (operating on angles specified in radians) and 2○ right argument gives the cosine function.
○1 (pi)
3.14159
○2 (2×pi)
6.28319
1○○.5 (sine of pi÷2)
1
2○○.5 (cosine of pi÷2)
6.13658E¯17
'1○X←' SMP_INT 0,(○1),100 (integral of sine from 0 to pi)
2
For a final example, consider calculating areas under the Standard Normal probability density function. First to define NORMDEN:
Hence, to calculate the area under the curve from 0 to 1:
and so the area under the curve corresponding to mean plus or minus one standard deviation is approximately 68%. Similarly,
thus the area to the right of 1.96 is .5-0.475, or .025
Producing Normal Distribution Tables¶
With a little program to loop around an array of appropriate values, tables of the Standard Normal distribution can be constructed easily, though it may take some time, depending on the speed of your computer. Here is such an example:
∇ R←DIM NORMTAB G;C
[1] R←(⍴G)⍴0 ⋄ C←1
[2] lp:R[C]←0.5-'NORMDEN' SMP_INT 0,G[C],100
[3] →((⍴G)≥C←C+1)/lp
[4] R←DIM⍴R
∇
As an example of its use consider the rather coarse grid of values from 0 to 2.8 in steps of .2
R←⎕←3 5 NORMTAB ¯1↓GRID 0 3 15
0.5 0.4207 0.3446 0.2743 0.2119
0.1587 0.1151 0.08076 0.0548 0.03593
0.02275 0.0139 8.198E¯3 4.661E¯3 2.555E¯3
With the help of some labelling programs (see the separate article) a table suitable for printing can easily be produced.
'Z value/.0/.2/.4/.6/.8/' CLABEL '0/1/2/' RLABEL ⍕R
Z value .0 .2 .4 .6 .8
------- -- -- -- -- --
0 0.5 0.4207 0.3446 0.2743 0.2119
1 0.1587 0.1151 0.08076 0.0548 0.03593
2 0.02275 0.0139 8.198E¯3 4.661E¯3 2.555E¯3
These calculations yield results which agree with commercially available tables to four decimal places.
Monte Carlo Simulation for Integrals¶
Finally in this article, consider the possibility of Monte Carlo Simulation to evaluate integrals. This is a technique of sampling the function f to be integrated uniformly at random throughout the range of integration and calculating the resulting sample mean. This is an approximation to the expected value of the random variable f(X), (where X has a Uniform distribution on [a,b]) which in turn differs from the integral required by a factor (b-a).
This approach is easily done in APL with the aid of a few more tools. First of all, generate values uniformly at random through the interval a to b:
∇ R←N UNIFORM A_B;A;B;T
[1] A←A_B[1] ⋄ B←A_B[2] ⋄ T←10000000000
[2] R←A+(B-A)×(0.5+?N⍴T)÷T
∇
10 UNIFORM 0 2
1.9156579 1.192201 0.097855169 0.015331741 0.50314947
0.19924829 1.4752251 1.70138 1.8358242 0.75769372
Add to the toolkit a program to calculate the mean:
Now calculate twice the mean of a set of function values (twice because the interval is of length 2):
Again, such a procedure can be called from a program, called MC_INT
∇ ∆R←∆FN MC_INT ∆A_B_N;∆A;∆B;∆N
[1] ∆A←∆A_B_N[1] ⋄ ∆B←∆A_B_N[2] ⋄ ∆N←∆A_B_N[3]
[2] ∆R←(∆B-∆A)×MEAN⍎∆FN,' ',(⍕∆N),' UNIFORM ∆A_B_N[1 2]'
∇
'FUN' MC_INT 0 2 1000
0.58683034
Of course, repeating the same command will give different estimates of integral. Moreover if you play around with it you will soon see that you have to use a lot of values to get a decent accuracy. So the technique is not used for simple integration. It does however have its uses in more complex situations.