Essays/Newton's Method
Newton's method is used to find approximations to a root of a differentiable function.
Introduction
N=: 1 : '- u % u d. 1' NB. prior to J9 require 'math/calculus' NB. J9 onwards N=: 1 : 'y-(u y)%0 u sslope_jcalculus_ 1]y' NB. J9 onwards (_2 + *:) N 1 1.5 (_2 + *:) N^:2 ]1 1.41667 (_2 + *:) N^:3 ]1 1.41422 (_2 + *:) N^:_ ]1 1.41421 2 - *: (_2 + *:) N^:_ ]1 4.44089e_16
N is an adverb where u N is one iteration for finding a root of u , whence u N^:n x is the result of n iterations on an initial estimate of x and u N^:_ x is the limit (to within the comparison tolerance) of the iterations.
Convergents
The convergents (results of iterations) obtain with an appropriate argument to the power operator ^: .
(_2 + *:) N^:(i.8) 1 1 1.5 1.41667 1.41422 1.41421 1.41421 1.41421 1.41421 (_2 + *:) N^:a: 1 1 1.5 1.41667 1.41422 1.41421 1.41421 2 - *: (_2 + *:) N^:(i.8) 1 1 _0.25 _0.00694444 _6.0073e_6 _4.51061e_12 _4.44089e_16 4.44089e_16 _4.44089e_16
Rational Numbers
If u uses only rational operations, then the iteration may be modified to produce rational results for a rational initial estimate. In such cases use of _ or a: in ^: should be avoided as the function may not have a rational limit.
To achieve this rational version of "N", we will create local copies of two of the calculus functions we loaded with "require 'math/calculus'". We keep sslope the same:
sslope=: 2 : 0 assert. _1 = 4!:0 <'m' NB. u is verb assert. 0 = 4!:0 <'n' NB. v is a noun assert. 0 = #$n NB. n is an atom assert. n > 0 NB. n > 0 (u. derivsecant n)"(u. f. b. 0) )
However, we need to modify its sub-function derivsecant by allowing for replacing the floating point epsilon 1e_7 with a rational equivalent 1r10000000:
derivsecant=: 2 : 0 : NB. Get function to use: u. or higher-order secant if. n = 1 do. func =. u.@] else. func =. u. derivsecant_jcalculus_ (n-1) end. NB. x must be an atom or conform to shape of y if. 0=#@$x do. x =. ($y)$x end. NB. replicate atom assert. x -:&$ y NB. shapes must agree
x =. x + ((1e_7"0)`(1r10000000"0) @. (64 128 e.~ 3!:0) y) * 0 = x NB. replace 0 by epsilon; preserve rationality.
newy =. y +"(#@$y) (,~$x) $ (#~ 1 j. #) ,x NB. array of moved points, each an array with 1 nonzero f0 =. x func y NB. the function at the initial point ((x func"(#@$y) newy) -"(#@$f0) f0) % x NB. evaluate function at moved points, calc slope. x used only for higher orders)
This bases the rationality of the epsilon on the rationality of y. This lets us define our rational version of Newton's method using the new local copy of sslope like this:
Nr=: 1 : 0 1r10000000 u Nr y : y-(u y)%x u sslope 1]y )
This allows us to control whether or not we get a rational result by applying the function to a rational argument, thus allowing us to extend the precision of our result arbitrarily:
(_2 + *:) Nr^:(i.8) 1x 1 30000001r20000001 17000002000000080000001r12000001400000060000001 ... 0j50 ": ,. (_2 + *:) Nr^:(i.8) 1x 1.00000000000000000000000000000000000000000000000000 1.49999997500000124999993750000312499984375000781250 1.41666666805555497685187770061704603904138803679948 1.41421568636341796166017977474869799833814519845978 1.41421356237476513837948736696017201172977125094212 1.41421356237309504886073629153533495244833341821383 1.41421356237309504880168872629734476818780034087328 1.41421356237309504880168872420969815237912581334999
Here we calculate the square root of 2 to an arbitrary precision by solving . We can see how much closer we get with each iteration:
2 - *: (_2 + *:) Nr^:(i.8) 1x 1 _99999979999999r400000040000001 _1000000799999919999969999997599999919999999r144000033600003400000192000006400000120000001 ... 0j_5 ": 2 - *: (_2 + *:) Nr^:(i.8) 1x 1.00000e0 _2.50000e_1 _6.94445e_3 _6.00756e_6 _4.72373e_12 _1.67012e_19 _5.90476e_27 _2.08765e_34
Here we see that seven iterations gives us about 33 digits of precision.
See also
Contributed by Roger Hui. An earlier version of the text appeared in the J Forum on 2007-10-13.