Fifty Shades of J/Chapter 23
Table of Contents ... Glossary ... Previous Chapter ... Next Chapter
Principal Topics
- | (magnitude / residue), |. (reverse / rotate), D. (derivative), {. (head / take), ] (same / right), ^: (power conjunction), _ (infinity), { (from), {: (tail), }. (behead / drop), -: (halve), #. (base), >: (increment), inner product, recursion
Introduction
In this essay it is assumed that the reader has at least a basic interest in, or knowledge of, elementary numerical analysis, and in particular of the three requirements of linear interpolation, root-finding, and numerical integration.
Linear Interpolation
Suppose a curve is known to go through the points and . The point at which the straight line connecting these points cuts the x-axis is , where . This computation is the basis of the simplest technique for finding roots. In the verb lint defined below, the left argument is x0 x1, while y0 y1 is the right argument.
lint =. dyad : '((|. x) -/ .* y) % -/ y' 1 3 lint _1 2 1.66667
Root-finding : Newton’s Method
Three functions will be used as as test cases throughout the rest of this article. These are
f1 =: #.& 1 _1 _2 NB. polynomial x^2 - x - 2 NB. base is ingenious, however f1 =: _2 _1 1&p. NB. is straightforward. NB. if f1 must use base, it should be NB. f1 =: (#.&1 _1 _2)"0 NB. otherwise f1 3 6 7 8 results in a confusing length error, NB. and f1 1 2 3 just gives something weird. f2 =: 2&o. - (*^) NB. cos(x) – x * exp(x) f3 =: 2 0 0 _1&p. NB. 2 – x^3
These test verbs will appear as a left argument to the numerical method adverbs. First Newton’s iterative formula for which the right argument is the initial guess. Newton is an adverb qualifying the function whose root is to be found:
Newton =: adverb : ']-u%(u D.1)'(^:_)("0) NB. pre J9
require 'math/calculus' NB. J9 onwards
New1=: adverb : 'y-(u y)%0 u sslope_jcalculus_ 1]y' NB. J9 onwards
Newton=: adverb : 'u New1 ^: _ ] y' NB. Iterate to convergence
f1 Newton 1
2
f2 Newton 1
0.517757
f3 Newton 1
1.25992
Root-finding : Secant Method
This requires two initial guesses between which the root is to be found. The result is the two most recent approximations. Like Newton, it is an adverb qualifying a function. 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 values to kick off the next approximation. This means that the infinity (convergence) option with the power conjunction is not available, instead the number of iterations is given explicitly as a left argument.
sec =: adverb : '}.@] , ] lint (u every@])' (^:[) 7 f1 sec 0.5 1 1.99963 2 4 f2 sec 0.5 1 0.517767 0.517757 6 f3 sec 0.5 1 1.25991 1.25992
To clear up the untidiness of receiving two result values when only one is really wanted, here is a short iteration verb which stops when things get sufficiently close :
secant =: adverb define
:
sec =. }.@] , ] lint u every@]
t =. sec y
while. (x < | t - r =. sec t) do.
t =. r
end.
{. r
)
0.00001 f1 secant 0.5 1
2
0.00001 f3 secant 0.5 1
1.25992
Root-finding : Regula Falsi Method
This is a variation of the secant method, sometimes called the Method of False Position, in which the end points at each step are reordered to ensure that the function values at the next iteration are of opposite sign. The left argument is the precision at which iteration is to stop. Ken Iverson described this in Math for the Layman
regf =: adverb define
:
z =. u every y
i =. y lint z
if. +./ x > | i - y do.
return.
end.
yi =. u i
if. yi =&* 1 { z do.
r =. x u regf (0 { y) , i
else.
r =. x u regf i , 1 { y
end.
)
0.00001 f1 regf 0.5 1
1.99999
0.00001 f2 regf 0.5 1
0.517755
0.00001 f3 regf 0.5 1
1.25992
http://jsoftware.com/pipermail/programming/2018-September/051922.htTacit adverb for false position methodml
Root-finding : Illinois Method
This is a variation adaptation of the previous method which was described in Vector 12.2 pp. 87-94 and Vector 12.4 pp. 98-106. It compensates for a curve turning more steeply on one side of an interval by halving the relevant ordinate. A second iteration step is inserted into regf :
illi =: adverb define
:
z =. u every y
i =. y lint z
if. +./ x > | i - y do.
return.
end.
yi =. u i
z =. z * 0.5 1 |.~ yi =&* 0 { z
i =. y lint z
yi =. u i
if. yi =&* 1 { z do.
x u illi (0 { y) , i
else.
x u illi i , 1 { y
end.
)
0.00001 f1 illi 0.5 1
_0.999971
0.00001 f2 illi 0.5 1
0.517755
0.00001 f3 illi 0.5 1
1.25988
There is no particular magic about the 0.5, and 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 a result where simpler methods fail. There is no universal best root finding function, and the choice of method in any particular case belongs to the realm of the Numerical Analyst’s art.
Integration : Simpson’s Method
The following methods require the function to be assigned to fn in advance of the call to the verb. The right argument is the range of integration. Simpson’s method is the basic way of doing numerical integration when analytic methods fail. The left argument is the number of intervals into which the range must be divided, and must be an even number. The greater the number the greater is the accuracy of the approximation to the integral.
Simpson =: adverb define
:
x =. >.&.-: x NB. ensure x is even
h =. - y -/ . * % x
(h*+/(1,((x-1)$4 2),1)*u every ({.y)+h*i.x+1)%3
)
2 f2 Simpson 0 0.5
0.303737
4 f2 Simpson 0 0.5
0.303783
4 f1 Simpson 0 0.5
_1.08333
4 f3 Simpson 0 0.5
0.984375
Integration : Adaptive Method
This is a development of the Simpson method which subdivides the range of integration in two at each step, and carries on applying Simpson integration recursively in each half until the required level of approximation specified in the left argument is attained.
adapt =: adverb define
:
r =. 4 u Simpson y
if. x > | r - 2 u Simpson y do.
return.
else.
m =. 0 1 |. each y , each -: +/ y
+/ (x u adapt {.&> m) , x u adapt {:&>m
end.
)
0.00001 f3 adapt 0 0.5
0.984375
0.00001 f1 adapt 0 0.5
_1.08333
0.00001 f2 adapt 0 0.5
0.303786
Integration : Romberg’s Method
The next variation is the Romberg method which uses series expansions for the error terms in Simpson approximation 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 =. adverb define
:
r =. > (2 ^ 1 + i. x) u Simpson every < y
while. 1 < # r do.
r =. }. r + (_1 + 2 ^ 4 * 2 + x - # r) %~ r - _1 |. r
end.
)
Romberg =: adverb define
:
t =. (i =. 1) u romb y
while. x < | t - r =. (i =. >: i) u romb y do.
t =. r
end.
)
0.00001 f2 Romberg 0 0.5
0.303783
0.00001 f1 Romberg 0 0.5
_1.08333
0.00001 f3 Romberg 0 0.5
0.984375
These few verbs and adverbs make root-finding and integration both practical and self-describing, and thereby give some insight into the routines offered within the more sophisticated mathematical packages.
Code Summary
Linear Interpolation
lint =: 4 : '((|. x) -/ .* y) % -/ y'
f1=: _2 _1 1&p. NB. polynomial x2-x-1 test case 1 f2=: 2&o.-(*^) NB. cos x – xex test case 2 f3=: 2 0 0 _1&p. NB. 2 – x3 test case 3
Root-finding
Newton=: 1 : '- u % u d. 1' NB. prior to J9 require 'math/calculus' NB. J9 onwards Newton=: 1 : '- u % u deriv_jcalculus_ 1' NB. J9 onwards
Newton=: (1 : ']-u%(u D.1)'(^:_))("0) NB. prior to J9
New1=: 1 : 'y-(u y)%0 u sslope_jcalculus_ 1]y' NB. J9 onwards
Newton=: 1 : 'u New1 ^: _ ] y' NB. Iterate to convergence
secant =: 1 : 0
:
sec =. }.@] , ] lint u every@]
t =. sec y
while. (x < | t - r =. sec t) do.
t =. r
end.
{. r
)
regf=: 1 : 0
:
z =. u every y
i =. y lint z
if. +./ x > | i - y do.
return.
end.
yi =. u i
if. yi =&* 1 { z do.
r =. x u regf (0 { y) , i
else.
r =. x u regf i , 1 { y
end.
)
illi=: 1 : 0
:
z =. u every y
i =. y lint z
if. +./ x > | i - y do.
return.
end.
yi =. u i
z =. z * 0.5 1 |.~ yi =&* 0 { z
i =. y lint z
yi =. u i
if. yi =&* 1 { z do.
x u illi (0 { y) , i
else.
x u illi i , 1 { y
end.
)
Integration
Simpson=: 1 : 0
:
x =. >.&.-: x NB. ensure x is even
h =. - y -/ . * % x
(h*+/(1,((x-1)$4 2),1)*u every ({.y)+h*i.x+1)%3
)
adapt=: 1 : 0
:
r =. 4 u Simpson y
if. x > | r - 2 u Simpson y do.
return.
else.
m =. 0 1 |. each ([echo)y , each -: +/ y
+/ (x u adapt {.&> m) , x u adapt {:&>m
end.
)
Romberg =: adverb define
:
t =. (i =. 1) u romb y
while. x < | t - r =. (i =. >: i) u romb y do.
t =. r
end.
)
romb =. adverb define
:
r =. > (2 ^ 1 + i. x) u Simpson every < y
while. 1 < # r do.
r =. }. r + (_1 + 2 ^ 4 * 2 + x - # r) %~ r - _1 |. r
end.
)