Clojure Dojo 3: From Heron to Newton-Raphson

;; We've so far written a square root solver using Heron's method

(defn average [a b] 
  (/ (+ a b) 2))

(defn abs[x] 
  (if (< x 0) (- x) x))

(defn make-good-enough? [n]
  (fn [guess] (< (abs (- n (* guess guess))) 1e-6)))

(defn make-improver [n]
  (fn [guess] (average guess (/ n guess))))

(defn iterative-improve [x improve good?]
  (if (good? x) x
      (iterative-improve (improve x) improve good?)))

(defn square-root [n]
     (iterative-improve 1.0 (make-improver n) (make-good-enough? n)))



;; Now in fact, Heron's method turns out to be only a simple example of a more general method of root finding known as Newton-Raphson.
;; Their great idea was to use the derivative of a function to help find its roots.

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Finding Derivatives (by cheating)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; Suppose we've got a function, say cube

(defn cube [x] (* x x x))


;; And we want to find out what its derivative is, i.e. how much it changes when you change the argument.

(cube 4)        ; i s 64
(cube 4.000001) ; is 64.000048000012

;; In other words, if we add 0.000001 to x, then (cube x) goes up by 0.00004800..... which is about 48 times as much.

;; We say that the derivative of cube at 4 is (about) 48.

;; It's different in different places.


(cube 3)        ; is 27
(cube 3.000001) ; is 27.000027000009005

;; So the rate of change of cube at 3 is (about) 27.

;;In general, we want a function that takes a function, and gives back a function that tells us how much it changes if you add a tiny bit to its argument.

(defn deriv [f]
  (fn [x] (/ 
           (- 
            (f (+ x 0.000001)) 
            (f x))
           0.000001)))


;; Let's try that out:

((deriv cube) 4)       48.00001200067072
((deriv cube) 3)       27.00000900546229

;; Notice what we did here! This is a function which takes a function and gives back a function.

;; We could name the answer!
(def dcube (deriv cube))

(map dcube (range 10))        ;(1.0E-12 3.0000029997978572 12.000006002210739 27.00000900546229 48.00001200067072 75.00001501625775 108.00001800248538 147.00002100198617 192.00002384422987 243.0000268986987)


;; So now we know how to calculate (a good approximation to) the derivative of any function.


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Newton Raphson
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; Imagine that we wanted to find the cube root of 100.

;; That's the same as finding a number which, when you cube it and subtract 100, gives you zero


;; So we should take as our equation to solve: (cube x) - 100 = 0

;; or (f x) = 0 where f is 

(defn f[x] (- (* x x x) 100))


;; Isaac Newton told us that if you have an equation like that, and you have a guess at a number
;; which will make it zero, then you can find a much better guess if you know the derivative.


(def df (deriv f))

;; Newton said: Suppose that I guess that 4 is the cube root of 100.

;; Then I try it (f 4) = -36

;; That's not so close. What's the derivative there?

(df 4) 48.00001200067072

;; What that means is that when we raise x by a tiny amount, say 4.0000001, then (f x) goes up by about 48 times that.


;; So Newton tells us:

;; The function is too low by 36.
;; If we make the guess a bit larger, then the function will go up by ~ 48 times as much.
;; We should try adding 36/48 as our next guess.

;; Let's try

(+ 4.0 (/ 36 48)) ;4.75

(f 4.75)          ;7.171875

;; ok, better. We were under by 36, now we're over by 7. 4.75 is too high!

;; Take the derivative again:

(df 4.75) ;67.68751426022845

;;Divide the amount we're off by the derivative
(/ (f 4.75) (df 4.75)) ;0.10595565634789487


;; So our next guess should be (- 4.75 0.10595565634789487) , which is 4.644044343652105

;; Let's try 4.644 :

(f 4.644) ;0.15592198400001678 ...homing in...

(- 4.644 (/ (f 4.644) (df 4.644)))  ;4.641590085793267

(f 4.64159) ;7.538717167676623E-5 ...very good...


;; Obviously we've got another guess-improving function here

(defn improve-cbrt100 [guess]
  (- guess (/ (f guess) (df guess))))


(improve-cbrt100 4.0)     4.749999812489567
(improve-cbrt100 (improve-cbrt100 4.0)) 4.644044335287887

(take 10 (iterate improve-cbrt100 4.0))   ; (4.0 4.749999812489567 4.644044335287887 4.6415901322397985 4.641588833613422 4.641588833612778 4.641588833612778 4.641588833612778 4.641588833612778 4.641588833612778)


;; Let's try 4.641588833612778, which seems to be about as good as floating point arithmetic can get us.

(f 4.641588833612778)      ; -2.8421709430404007E-14

(cube 4.641588833612778)   ; 99.99999999999997

;; Pretty good for only 5 steps!

;; What about a good-enough function to tell us when to stop iterating?


(defn good-enough-cbrt100? [x]
  (< (abs (f x)) 0.0000001))


;; Now we can plug the improver and the good enough function into iterative-improve, as before. We'll use 1.0 as our initial guess again.
(iterative-improve 1.0 improve-cbrt100 good-enough-cbrt100?) ;; 4.641588833613406

(cube 4.641588833613406)   ;; 100.00000000004056



;; But of course, nothing we did depended on (f x) being (- (cube x) 100)
;; We can make a guess-improver for any function



(defn make-improver [f]
  (fn [guess] (- guess (/ (f guess) ((deriv f) guess)))))


;; and a function to tell us whether it's good enough

(defn make-good-enough [f]
  (fn [guess] (< (abs (f guess)) 0.0000001)))

;; What if we'd like to solve the equation x^3 + x^2 + x + 1 = 0 ?


(defn solve [f guess]
  (iterative-improve guess (make-improver f) (make-good-enough f)))

(solve (fn [x] (+ (* x x x) (* x x) x 1)) 1.0) ;  -1.0000000235152005

;; Let's see how good an answer that is.

((fn [x] (+ (* x x x) (* x x) x 1)) -1.0000000235152005) ; -4.703040201725628E-8




;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Newton Raphson Solver
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; Here's the whole of our general function solver, but let's make the two 'magic numbers' into variables too.

(defn deriv [f dx]
  (fn [x] (/  (- (f (+ x dx)) (f x) )  dx)))

(defn make-improver [f dx]
  (fn [guess] (- guess (/ (f guess) ((deriv f dx) guess)))))

(defn make-good-enough [f tolerance]
  (fn [guess] (< (abs (f guess)) tolerance)))

(defn iterative-improve [x improve good?]
  (if (good? x) x
      (iterative-improve (improve x) improve good?)))

(defn solve [f guess dx tolerance]
  (iterative-improve guess (make-improver f dx) (make-good-enough f tolerance)))


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Some applications
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;usually we'll want an initial guess of 1.0, a tolerance of 0.0000001 and a dx of 0.0000001
(defn solve-it [f] 
  (solve f 1.0 0.000001 0.000001))

;; Square roots are the solutions of (- (* x x) n) = 0


(defn sqrt [n] 
  (solve-it 
   (fn [x] (- (* x x) n))))

(sqrt 2) 

;; Cube roots are the solutions of (- (* x x x) n) = 0

(defn cubert [n] 
  (solve-it 
   (fn[x] (- (* x x x) n))))

(cubert 2)


;; General positive integer powers can be expressed recursively:

(defn pow [x m]
  (if (= m 0) 1
      (* x (pow x (- m 1)))))

(pow 2 3) ;8

;; And if we can express them, we can use our method to find mth roots.


(defn mthrt [m n] 
  (solve-it 
   (fn[x] (- (pow x m) n))))

(mthrt 4 10000) ; 10.0
(mthrt 123 234) ; 1.045350466874579

(pow 1.045350466874579 123) ;; 234.0000000243792


;; But of course, we can solve any equation which is amenable to Newton's method.
;; What number, when multiplied by itself 5 times, and added to its cube, is equal to 700?

(solve-it (fn[x] (+ (* x x x x x) (* x x x) (- 700))))         ; 3.653807138770812

((fn[x] (+ (* x x x x x) (* x x x))) 3.653807138770812) ; 700.0000003041337

;; And if we'd like a better answer, we can tighten the tolerances

(solve (fn[x] (+ (* x x x x x) (* x x x) (- 700))) 1.0 0.0001 0.0000000001) ; 3.6538071384442974


((fn[x] (+ (* x x x x x) (* x x x))) 3.6538071384442974) ;  700.0000000000817









;; In theory, we'd be able to use exact arithmetic to get arbitrarily accurate answers.
;; Unfortunately, there's a problem:


(def it (make-improver (fn[x] (+ (* x x x x x) (* x x x) (- 700))) 1/1000000)) ; let's use an exact number for the derivative approximation

;; and start with an exact guess

(it 1)                  ; 706000013000011000005000001/8000013000011000005000001
(it (it 1))             ; 701635321996143769309418549344377227910245191235491382603948489939699061403896688219959595828904782223381164841011935230760200930031084004228000441000030000001/9938320550977024134920249980653903204519665284470123223646489755187942893061019955878983167997749088761731121640594804980166030029688004228000441000030000001

(it (it (it 1)))        ; 680237917441885524336418138840317390387477847859858876595840234272149425820646224687954143292874376454521646334951672503154114880181300843385953617515081285781202235479686701822851428206060076054501437170597117295185786622448092483300621863099620331252416621986363001610323019591307968502301991515214451683404838437583955509995832608492690263029459069929976980341191436712932448855080982518511841862770021155980758185192008036141486759432138585947862423162360655529584841649322416573707879742014567586254212003318865179086757633089843219756812844011861068141391297639719428726706532258163565342026653477095733524189143027946838723155848407864793520497990564324241350635343653518356286709120797176348156180962741126448716721627967605789740379911361208362221973935469409915879343921377079168493057613428011966000155000001/12044286332140129752156903683583126447244938902908926261607579552042467607850219073475448382767207875539137993119265514088232278316554495937166615076210884321449749330293196083213705419806298445939325718843197812663753008605821508903073125871316024961234609794639656768170967982087039663880189123415565181043003755509161080647377737131911463013608843269269085395872146930064532858794789612245264182978492519500766787287053880913766446555513734299742826803443919934807497308923876636662216418140448109767596616697192628702623509472332595967524822599329261416093615390398256303064181652715052013431614546617380753560488219680455031256525883523608904103537508019124887595921116604264464741508409082509232101707795378534258208496424946224563887570504744604783840489200142804159760458106765068490963613428011966000155000001
(it (it (it (it 1))))   ; 582683064888886254627237413128914355133205070107730345306566818667169839077043977063385150567627868898954627396888803873298542915735291068231854737660749871284970969691761755237023440170586132711096078205774661890184691373273510273283202351760392701768411497903999874829223861118399238720882465581354352630000971969184626801423877016182068328118391265021597001152296290812122250720956031547021906861678543752400823983366892551426008467510890046723277161664312028285174276774217375277065253742530859064530016440400316790801769530544371284081707010989839631883556283241113998175234840943170722268328720701199347867588322060049770935643979286527454236006759544836737208537749694757872416290454093930058675287658912516799102545283953465977040744182800487171767210596289815719837983663631950117171893161135122746529146479448859227394441049923190055342857801930221659120398700765917475396756356933748228771216329690971598121024906393909199740001655437226675292089740865401850681717862593295044249444227096530596676969537355119944076825364494231692463504463385356777816358752007019352340833078969061391374565063342370784189339508167803312120387247644380450529051088582048118902223778854166538058459895429188150443034952695968773967230092326545545925098065274390529306120532407597280407465373342806740914557404279803404281875098412568052206856189504731510965158306039056809060766012400914292618679128113015139253847973084750652345634913827827043851034605686643123187076375915598639958363937697829956243245104853978610467265918061529051311071288413580756443693449482971629647374397753310557728689142976074021187362385338932147567642210185904671961035565315840293902065425141448705637965586698203041117408261444705945865047531954486514780505078832025358426644461672866572253641788067682868313500867194340509941792198538459647737009443599755837217723308620719398896201880698984796565219264280137596609642845400547125754456209756196169964659473873586279764214453671645956003619308063592648337612774637171241884690676604002029522916589191903605327335495103912144305652194435971484948567480613108531306049822517644426959900990535402937677569195079220088350105434599994075059763835418492356744894787315782423704801923419774013400657508864679447262609755156219146757585536473633170206959340454211578712959760853352548696411404380901268846622462357056396277653709161131992015278410666198586576834004214031859566710336893361077448333057320403538099709665120727817043295276197596679826722665947871368777705493013687853852488123574246109063974638016285172124627788918704893034469219280568313768497634816004319733836637800496141268248807282891107795078010343054429751145414810367778479411850027960387501665140659520089944949565082934980430105203072555206594274562501052867096198848952152330124216130716659045181152784334507386057959224529758674635647399911089583876402241119812642070687710563501331168402312038757841360479166547376873053789334438351535437225196907596309808445013472858435271020713482531473249425954713138213993853690048279402486411516169465400907740591254918376665681348518024822387460412485718766760961958943780464278973812982576607564820484862661878406223946061213224303753767486290224002094709777464746472645414893385625996401380435392614667079096134453611792569645198318837308794387253480934803497220222891193495355138449277282306273087395040974247699578508944806995861199096983093778899088019827135297023837526798957092191839173362464128894165682017815738236599955940542605196688759536725013516191327597134380086354086023426735965448551045630879574037331690144148829216502025494635660367861837832650027912461998002111883455796634997105423180569440930861025594379890711775754266620157029760796036618547810092105837604672377921726136778076482261042587195382143724289727648562704136017853509498366890447377518135820019350488687720760940855188770925017427141691614329578171795155245502836823003566045473401512439769020488822239922252814521793930372045109677392922267054560823873124298362347797278136358895228352710628394677753104291145730759681035568127201415728926581298819651031969039124011747083909428303966000780000001/12896626453600069493698610007462574074360012803192517550652653321288458704437628147812351238096058570880646809642684891850066999395344051036048084092196565019966845694603226043392692831832408470853928980193581747252243200468103361176860667342346033755639391767792595834026685821301417812313001324714923798171069125692865759809093740132786844460822717151100184874113026243342113737733017773505307764700961069595789996009548181632271209787338405203417470570243624108296656484641803537134016638223204512186999071878554046958682488174902216383193717986143532189831044292308519842126873269313853560876247133978458353574480851312759631684358889032900084642745496288725888743073969677137896558829603592262230802522695739419744389018487062275548598153935413565238867472051597918163575900094811377931457308002870398147828447466466856465628782660404829455110767346686507193439951612181709790999494520905843712219628937059390531846716829881905679762662797797772057202910931972704357372155619580235092749976889122762716133368842037661984150683629394903567236216963194674605967568392431208794043431169892910445094675188563758743114907987788746591274250207082015901133815715910101081963885658315948409184499380715077775262542275170959223244947248063594716133379982128239636514206935394270725149984811164184816315321673225284930879745768846720157626335366782084507443322775922411153536241445645457909068808146547276992161599402164331016798612763537783551540734336332823264126557187267639802943507420912300582217231591699364446110633032192667875702437386209756724502546707178905247996448449229725676664533990896268840252253987954869584921155755185596659051602784734099723461617865205317857691298850939510647703168248071840078142301962876095991821410664592847145496334778994271130316719075179553320153706951867058953414639227099786449755622188371176105537376510675813747465892564129310902573333355338988514058192477421861710533925137494638802746649014987114280003213204406857866359359280991867594433240560958642355778026603456310571266226381230212501928363154073509767737834411081511741666078389988426255700754713921887266615912226559289922110288321637513352794999588449807206908304314944589521175612052878625277657474456489493344970961101196034555416517832906272790448598835563702334203053761460216483188413922273789078424479313993102344574019369085117678979346999182040112416824764512591041037973913657403354399232242224999268943921208792520828958494399634230250208344108593016608457066759423447434130404474515620954899757483291938447043440289286090559545316946393298710577732756942628616513884579585056022530728783711122920500882377115985715201353863643271746094708630198929184352363067785421438134909914864758356500163689964428638647687495591442052889991128729671129301459039511338151733777970655351650119214267552140948065933092934877227143743916767768243644993862805442252649852133799824872693054837337674268061282267623038576543034125243858631285899198303109966963172209152874196720199801119238717086988413507996967711063340008844027377086930137204987454203538469697377749646520851213248589905149021586496254737543585289159013650128559023578482611563124559893592352427180734048684239663086639329223310790701294430611951452851943347414551549066596198403402509218874978312155181397196069302152786264584396523101878569322693955950044160266238264131771759791930583114290189937813208245285469977694328792320835576523730805116731577483465467653183154605688937756438221438317995677209102959548247744261605934633001969642392983652820479907300076279299375999438660224810860948463249531812666274785951673129729793010427319485797202888099118900435857053255162204718549364268309991788374329365868674979202327345545939405245168318555414785299112969636869024689405450765132794123797478424267092928772922594681937404357022580366085874534312247818159047202717252632019175403527618827567476810823460590807535175349230955481084676292524894761199878853286342204470938790343375932066863536221551379561495530913704293344747209001634639584668143339854848154541379814946625988808902291845280350306322621482114145676960211744291909428303966000780000001

;; after 4 iterations we've got 8000 digits in our fraction! It might take a while to calculate each step.

;; It's possible to control this crazed exponentiation of digits. 
;; I leave the problem of modifying the method to produce arbitrarily accurate answers as a problem to the reader.