RealInterval Newton Sample

This sample might look a bit more confused than usual as I will explain my thoughts while programming this example. I wanted to make a Newton zero-finder for intervals, since this is explained in every paper about interval arithmetic. First I thought, why should I look at these papers, everybody knows how this works and in the AutomaticDifferentiation documentation I did such a thing with five short lines. Unfortunately my quickly thrown together algo could not find any zero at all. After noticing my error, I thought following a simple how-to might perhaps not be completely wrong, and I then chose the algo decription in http://eprints.ma.man.ac.uk/1204/01/narep416.pdf page 27 which looks uncomplicated, and obviously invites further simplification. That was not a bad idea (you will additionally need AutomaticDifferentiation and DHB-wk for this):

precision :=1e-8.
(fp:=Fixpoint new) equalityTest: [:a :b| a < precision].
"results stores Arrays of the form: {aDualNumberOfaRealInterval. itsWidth} :"
results := Heap sortBlock: [ :a :b|
	(a at:2) >= (b at:2)] .	
fp block: [ :dummy|
	array := results removeFirst.
	array second < precision ifTrue: [results add:array] ifFalse:[
	x := array first. 
	n := x value midPoint -
             ((function value: x value midPoint)/(function value: x) eps).
	n := n intersection: x value.
	((function value: n) includes: 0) ifTrue: 
	 [n isIntervalUnion ifTrue: 
	   [x :=n first.
	    n :=n second.
	    results add: (Array with: (DualNumber value: x eps: 1) with: x width).
	    results add: (Array with: (DualNumber value: n eps: 1) with: n width).
	   ] ifFalse: 
	   [results add: (Array with: (DualNumber value: n eps: 1) with: n width)]
	 ]].
	array second] .

This is roughly the described algo in a short form that more or less works. The simplifications go a bit too far (eg perhaps the includes: 0 part). But then so what, it is better to start with a simple version that is easily understandable.

In AutomaticDifferentiation I had used this test-function:
f(x) = x³ + 2 x² + x
which has zeros at -1 and 0.

Let's test the algo:

results := Heap sortBlock: [ :a :b|
	(a at:2) >= (b at:2)] .
interval := -2 hull: 3.
results add: (Array with: (DualNumber value: interval eps: 1) with: interval width).
function := [:x|(x raisedToInteger: 3)+(2*x squared)+x].
fp evaluate.
results.-->"a Heap(
     an Array(DualNumber(value: [-0.9999999291210566,-0.999999925428417] eps: 1) 
              3.6926396385084104e-9) 
     an Array(DualNumber(value: [-1.2370295793295068e-9,1.5438236391280191e-9] eps: 1) 
              2.780853218457526e-9) 
     an Array(DualNumber(value: [-1.000000011933429,-1.000000011933429] eps: 1) 
              0.0) 
     an Array(DualNumber(value: [-1.0000002618674235,-1.0000002613573151] eps: 1) 
              5.101083999647926e-10))" 
results collect:[:a|a first value asNumber].
 "an Array([-0.9999999291210566,-0.999999925428417] 
     [-1.2370295793295068e-9,1.5438236391280191e-9] 
     -1.000000011933429 
     [-1.0000002618674235,-1.0000002613573151])"
results collect:[:a|function value: a first value asNumber].
 "an Array([-1.477056255083653e-8,1.4770551892695494e-8] 
           [-1.2370295793295068e-9,1.543823643894802e-9] 
           0.0 
           [-2.040502655731302e-9,2.0403658762546684e-9])"
results collect:[:a|function value: a first value midPoint].
 "#(-5.218048215738236e-15 1.5339702994631744e-10 0.0 -6.838973831690964e-14)"
r:=results collect:[:a| a first value midPoint].
 "#(-0.9999999272747369 1.5339702989925613e-10 
    -1.000000011933429 -1.0000002616123693)"

There are several noticeable things. First there are more results than one would want to have. Of course I could condense that to a more pleasing result with my usual, brash heuristic approach. Second, while 0 is indeed in the results, -1 actually has been thrown away. Third, the third result is a Number, not -1, but not too far away and f(-1.000000011933429) = 0.0 indeed. And that number in a way explains, why -1 has been lost:

function value:(DualNumber value: (r at:3) eps: 1)."-->
DualNumber(value: 0.0 eps: 2.3866858711585337e-8)"

In a way this is a perfect solution since f=0, it can't get any better. But its derivative is unequal to zero! And as mentioned in the AutomaticDifferentiation docu, the derivative calculation is rather exact. If this mini-interval would go just one more time through the iteration, it would get thrown away, because the next n would necessarily move outside the mini-interval with a width of zero.

In other words if you do heavy calculations (which produce floating point errors) and combine the results with set operations on small sets, you can get into trouble. Here one has more calculations than meets the eye, since the derivative calculation is done on intervals and every iterative solution deals necessarily with small sets (which multiply the importance of floating point errors).

I think this example shows nicely the shortcoming of this experimental kit: a real interval-arithmetic always encompasses the results irrespective of floating point errors, it takes those floating point errors into account and returns accordingly wider intervals. The other algos I described work better, because this is the first algo that uses interval-computation to numerically adjust a result (in a way move it via floating point errors) and then uses set operations on this result.

For another example let's find the zeroes of sin(x) for -4 ≤ x ≤ 9.

"here i had to lower the precision considerably (bigger precision number). 
if you raise the precision, you will loose solutions again"
precision :=1e-6.
results := Heap sortBlock: [ :a :b|
	(a at:2) >= (b at:2)] .
interval := -4 hull:9.
results add: (Array with: (DualNumber value: interval eps: 1) with: interval width).
function := [:x|x sin].
fp evaluate.
results. "a Heap(
     an Array(DualNumber(value: [3.1415926405844115,3.141592819059111] eps: 1) 
              1.7847469946374872e-7) 
     an Array(DualNumber(value: [6.2831853071583845,6.2831853082974] eps: 1) 
              1.1390151044565755e-9) 
     an Array(DualNumber(value: [-3.141592653616929,-3.1415926535897705] eps: 1) 
              2.715871971759043e-11) 
     an Array(DualNumber(value: [-6.606664114438068e-11,5.7965521488637095e-12] eps: 1) 
              7.186319329324439e-11))" 

results collect:[:a|function value: a first value].
 "an Array([-1.6546931768079595e-7,1.3005381782952002e-8] 
           [-2.1201951977932454e-11,1.117813152478643e-9] 
           [-2.277101033773543e-14,2.7135948707252694e-11] 
           [-6.606664114438068e-11,5.7965521488637095e-12])" 
results collect:[:a| a first value midPoint].
 "#(3.1415927298217614 6.2831853077278925 -3.14159265360335 -3.013504449775849e-11)"
results collect:[:a|function value: a first value midPoint].   
 "#(-7.623196817096688e-8 5.483060443395651e-10 
    1.355658884845748e-11 -3.013504449775849e-11)"

I guess that is enough of a demonstration that AutomaticDifferentiation and RealInterval can work nicely together, and that RealInterval's implemention, which does not take care of floating point errors, is occasionally rather dubious.

Back to Examples