J-ottings 19: Getting to the Root of the Matter
Vector 15:3, page 33
Transcribed from the printed issue; not yet reviewed. Read it in the PDF of the issue, page 33
Einstein is reputed to have said "God doesn’t care about our mathematical difficulties; he integrates empirically." The nonchalance with which thousands of scientists of all complexions, and even greater battalions of students, nowadays obtain numerical results involving huge amounts of computation prompts the observation that the combined skills of numerical analysis and computing science appear to have brought about at least some tiny encroachment on this aspect of the divine mind.
The new release of J contains a link to LAPACK (short for Linear Analysis Package) with which it is claimed that J joins the league of MATLAB, Mathematica, Maple etc. in providing access to such computing power, thereby countering the argument in Mathematics departments "Why should I concern myself with things like J when I have MATLAB?" This article concerns the many mathematically expressible problems which can be solved without recourse to the heavy artillery of MATLAB and LAPACK. Specifically, a root finder and an integrator can cope with a large range of tolerably well-behaved mathematical functions.
It is useful to start with a verb for linear interpolation between points \((x_0,y_0)\) and \((x_1,y_1)\). Its content should create little difficulty in understanding even for those with quite modest computer awareness:
linp=.3 : 0
:
r=.((|.x.)-/ .*y.)%-/y. NB. x. and y. are coords of 2 points
)
0 1 linp _0.1 0.5
0.1666667
An implicit, albeit less immediately comprehensible form, which is almost twice as fast is:
Now consider root finding. The differential operator D. allows the Newton-Raphson formula to be written as an adverb which looks very like its mathematical equivalent:
This describes just one step of the process, so use the power adverb to carry on until further iteration produces no change:
To test this and subsequent functions here are three simple verbs providing mathematical functions:
f1=.#.&1 _1 _2 NB. polynomial x2 -x -1
f2=.2&o.-(*^) NB. cos x - x.exp(x)
f3=.2&-@^&3 NB. verb to find the cube root of 2
f1 Newton 1
2
f2 Newton 1
_0.51774
f3 Newton 1
1.25992
The Secant Method is another root-finding technique. Since two points are needed to define a secant as opposed to one for a tangent, the secant adverb requires two start values, and at every step delivers a further two start values to define the next approximating secant. This means that the "infinity" option with power is not available, however it is still possible to give the number of iterations explicitly as a left argument:
sec=.1 : '}.@],] linp (x.&>@])' (^:[)
3 f1 sec 0.5 1 NB. 3 steps of Secant Method
1.4 1.66667
3 f2 sec 0.5 1
0.5158504 0.5177667
3 f3 sec 0.5 1
1.19838 1.24658
To clear up the untidiness of receiving two result items whereas eventually only one is really wanted, write a short iteration verb which stops when things get sufficiently close:
sec1=.1 : '}.@],] linp fn&>@]' NB. secant single step (monadic)
secant=.3 :0
:
t=.fn sec1 y.
while.(x.<|t-r=.fn sec1 t) do.t=.r end.{.r
)
eps=.0.00001
fn=.f1
eps secant 0.5 1
2
fn=.f3
eps secant 0.5 1
1.25992
A variation of the Secant Method, called Regula Falsi, or the Method of False Position, reorders the end points at each step in order to ensure that the function values at the next iteration are of opposite sign:
regf=.3 :0 NB. Regula Falsi=Method of False Posn
:
yi=.fn i=.y.linp y=.fn&>y. NB. i=interpolate, yi=fn. value at i
if. +./x.>|i-y. do.r=.i return. end NB. stop if successive i's close
if.yi=&*1{y do.r=.x.regf(0{y.),i NB. if yi,y1 of same sign choose y0
else.r=.x.regf i,1{y. end. NB. if yi,y0 of same sign choose y1
)
fn=.f3
eps regf 0.5 1
1.25992
As an aside on the J language, I find the verb =&* (equal with regard to sign) in the penultimate line is a pleasing application of the "with" conjunction.
A variation of the Regula Falsi Method called the Illinois Method was described in APL terms in Vector Volume 12 (see [1] and [2]). This compensates for a curve turning steeply at one side of the interval by halving the relevant ordinate. A second interpolation step is inserted into regf :
illi=.3 :0
: yi=.fn i=.y.linp y=.fn&>y. NB. start and finish as for regf
if. +./x.>|i-y. do.r=.i return. end.
if.yi=&*0{y do. NB. three lines inserted to include
yi=.fn i=.y.linp y=.1 0.5*y NB. an Illinois step which reduces
else.yi=.fn i=.y.linp y=.0.5 1*y end. NB. one of the ordinates by one half
if.yi=&*1{y do.r=.x.illi(0{y.),i
else.r=.x.illi i,1{y.end.
)
fn=.f1
eps illi 0.5 1
_0.9999968 NB. obtains different root from secant
fn=.f3
eps illi 0.5 1
1.25991
There is no particular magic about the 0.5, and presumably other fractions could be experimented with. The broad objective of such variations of the Secant Method is either to improve convergence for "difficult" functions, or sometimes to achieve it where simpler methods fail. There is no overall best root finding function, and the choice of method in any particular case belongs to the realm of the Numerical Analyst’s art, even if Mathematica, etc. is at his disposal. For well-behaved functions the variations given above rarely outperform the Newton or Secant methods in execution time.
Turn now to integration. As every schoolboy knows (provided he is in the Mathematics sixth!) Simpson’s rule is the basic way of doing numerical integration when analytic methods fail:
Simpson=.3 :0
:
(h*+/(1,((x.-1)$4 2),1)*fn&>({.y.)+(h=.-y.-/ .%x.)*i.x.+1)%3
)
fn=.f2
2 Simpson 0 0.5
0.3037372
4 Simpson 0 0.5
0.3037831
Not perhaps a J verb of compelling beauty, but the basic mathematics is reasonably recognisable. Two variations can be used to enhance it, first the adaptive method which recursively divides intervals in half and performs successive Simpson integrations in each half. If these agree within tolerance, stop the process for that interval, if not further divide the interval into halves, and perform the adaptive method in each half-interval.
adapt=.3 :0
:
r=.4 Simpson y. NB. Stop when first two Simpson
if. x.>|r-2 Simpson y. do. return. NB. estimates sufficiently close
else. m=.0 1|.each y.,each-:+/y. NB. Divide interval in two and
r=.+/(x.adapt {.&>m),x. adapt {:&>m NB. integrate in each sub-interval
end. NB. separately
)
fn=.f2
eps adapt 0 0.5
0.303786
The second variation is the Romberg method which uses series expansions for the error terms in Simpson estimation and applies corrections to progressively more refined Simpson estimates. The process is too complicated to describe in succinct comments, but can be found in most elementary texts on Numerical Analysis.
romb=.3 :0
: r=.>(2^1+i.x.)Simpson each<y.
while.1<#r
do. r=.}.r+(_1+2^4*2+x.-#r)%~r-_1|.r end.
)
Romberg=.3 :0
:
t=.(i=.1) romb y.
while.(x.<|t-r=.(i=.i+1)romb y.) do.t=.r end.
)
fn=.f2
eps Romberg 0 0.5
0.3037833
These few verbs and adverbs make root-finding and integration a practical proposition in many cases where the default approach might be to invoke MATLAB (if available!). Doing it this way is inherently satisfying, and offers insight into the beginnings of what lies inside the more sophisticated mathematical powerhouses.
References¶
- Dietrich Trenkler, Computing Clopper-Pearson Confidence Limits by the Illinois Method, Vector 12.2, pp. 87-94
- Norman Thomson, Technical Note on Confidence Limits, Vector 12.4, pp. 98-106