Arjen Markus (9 july 2003) Having read a book about numerical accuracy and stability, as well as several technical descriptions about interval arithmetic (certain Fortran compilers support it), I wanted to have this facility in Tcl too.
The idea: a number is not completely accurate, but there is some error margin. Then summing two such numbers becomes a matter of manipulating intervals:
x = [[a,b]] y = [[c,d]]
then:
x+y = [[a+c,b+d]]
for instance.
I know that it has severe drawbacks (and for these I want to provide some counter measures, like an alternative technique, called "running error analysis"):
Well, creating the procs was easy and fast, so we can focus on obtaining experience :) Here is the script with a simple example.
Note: the functional or prefix style of arithmetic is really useful here - it will allow us to overload the procedures without changing the code everywhere.
# interval.tcl --
# Rudimentary support for interval arithmetic in Tcl
#
# To do:
# Elementary functions - cos, sin, tan, abs, exp, log, pow
# General function - function
#
namespace eval ::math::interval {
namespace export + - * / \
contains interval begin-end \
square abs
}
# + --
# Add two interval numbers
#
# Arguments:
# operand1 First operand (interval number or regular number)
# operand2 Second operand (interval number or regular number)
# Result:
# Sum of the two intervals
#
proc ::math::interval::+ {operand1 operand2} {
foreach {begin1 end1} [concat $operand1 $operand1] {break}
foreach {begin2 end2} [concat $operand2 $operand2] {break}
return [list [expr {$begin1+$begin2}] [expr {$end1+$end2}]]
}
# - --
# Subtract two interval numbers
#
# Arguments:
# operand1 First operand (interval number or regular number)
# operand2 Second operand (interval number or regular number)
# Result:
# Difference of the two intervals
#
proc ::math::interval::- {operand1 operand2} {
foreach {begin1 end1} [concat $operand1 $operand1] {break}
foreach {begin2 end2} [concat $operand2 $operand2] {break}
return [list [expr {$begin1-$begin2}] [expr {$end1-$end2}]]
}
# * --
# Multiply two interval numbers
#
# Arguments:
# operand1 First operand (interval number or regular number)
# operand2 Second operand (interval number or regular number)
# Result:
# Difference of the two intervals
#
proc ::math::interval::* {operand1 operand2} {
foreach {begin1 end1} [concat $operand1 $operand1] {break}
foreach {begin2 end2} [concat $operand2 $operand2] {break}
set v1 [expr {$begin1*$begin2}]
set v2 [expr {$begin1*$end2 }]
set v3 [expr {$end1 *$begin2}]
set v4 [expr {$end1 *$end2 }]
set min $v1
set max $v1
foreach v [list $v2 $v3 $v4] {
set min [expr { $min<$v? $min : $v }]
set max [expr { $max>$v? $max : $v }]
}
return [list $min $max]
}
# / --
# Divide two interval numbers
#
# Arguments:
# operand1 First operand (interval number or regular number)
# operand2 Second operand (interval number or regular number)
# Result:
# Difference of the two intervals
# Note:
# The second interval may include 0!
#
proc ::math::interval::/ {operand1 operand2} {
foreach {begin1 end1} [concat $operand1 $operand1] {break}
foreach {begin2 end2} [concat $operand2 $operand2] {break}
if { $begin2 <= 0.0 && $end2 >= 0.0 } { error "Division by interval containing zero" } }
# contains --
# Check if the interval contains a given value
#
# Arguments:
# interval Interval number in question
# value Value that is or is not contained
# Result:
# 0 if not, 1 if it is contained
#
proc ::math::interval::contains {interval value} {
foreach {begin end} [concat $interval $interval] {break}
if { $begin <= $value && $end >= $value } {
return 1
} else {
return 0
}
}
# interval --
# Create an interval number from a centre and a tolerance
#
# Arguments:
# centre Central number for the interval
# tolerance Tolerance (relative error, defaults to 1.0e-10)
# Result:
# Interval number
# Note:
# For centre = 0.0, the tolerance is chosen to be absolute
# - a hack in fact.
#
proc ::math::interval::interval { centre {tolerance 1.0e-10} } {
if { $centre != 0.0 } {
set begin [expr {(1.0-$tolerance)*$centre}]
set end [expr {(1.0+$tolerance)*$centre}]
} else {
set begin [expr {-$tolerance}]
set end $tolerance
}
return [list $begin $end]
}
# begin-end --
# Create an interval number from two extremes
#
# Arguments:
# low Lower bound for the interval
# high Higher bound for the interval
# Result:
# Interval number
#
proc ::math::interval::begin-end { low high } {
if { $low < $high } {
return [list $low $high]
} else {
return [list $high $low]
}
}
# square --
# Return the interval containing the square of an interval
#
# Arguments:
# interval Interval to be squared
# Result:
# Encompassing interval
# Note:
# The trick is that the square of a number can not be negative!
# And the operands are correlated so that a direct multiplication
# gives a very pessimistic answer
#
proc ::math::interval::square { interval } {
foreach {begin end} [concat $interval $interval] {break}
set min2 [expr {$begin*$begin}]
set max2 [expr {$end*$end}]
if { $begin < 0.0 && $end > 0.0 } {
return [list 0.0 [expr {$min2>$max2? $min2 : $max2}]]
} else {
return [list [expr {$min2<$max2? $min2 : $max2}] \
[expr {$min2>$max2? $min2 : $max2}] ]
}
}
# abs --
# Return the interval containing the absolute value of an interval
#
# Arguments:
# interval Interval to be treated
# Result:
# Encompassing interval
#
proc ::math::interval::abs { interval } {
foreach {begin end} [concat $interval $interval] {break}
set min [expr {abs($begin)}]
set max [expr {abs($end)}]
if { $begin < 0.0 && $end > 0.0 } {
return [list 0.0 [expr {$min>$max? $min : $max}]]
} else {
return [list [expr {$min<$max? $min : $max}] \
[expr {$min>$max? $min : $max}] ]
}
}
# Test code
#
namespace import ::math::interval::*
puts [+ {0.9 1.1} {2.1 2.2}]
puts [+ {0.9 1.1} 1]
puts [+ 1 1]
puts [- {0.9 1.1} {2.1 2.2}]
puts [* {0.9 1.1} {2.1 2.2}]
puts [/ {0.9 1.1} {2.1 2.2}]
puts "Square:"
puts [square {1 2}]
puts [square {-1 2}]
puts [square {-1 -2}]
puts "Abs:"
puts [abs {1 2}]
puts [abs {-1 2}]
puts [abs {-1 -2}]
#
# Now a more practical example
#
puts "Integrating a parabolic function f(x) = alpha * x**2"
puts "Assumes that all steps have independent uncertainties!"
set alpha [interval 1 0.001]
set alphapl 1
set dx 0.1
set sum 0.0
set sumpl 0.0
for { set i 0 } { $i < 300 } { incr i } {
set x [* $i $dx]
set sum [+ $sum [* $alpha [* $x $x]]]
set sumpl [+ $sumpl [* $alphapl [* $x $x]]]
}
puts "(Interval) estimate of integral: $sum"
puts "(Plain) estimate of integral: $sumpl"The output:
3.0 3.3 1.9 2.1 2 2 -1.2 -1.1 1.89 2.42 0.409090909091 0.52380952381 Square: 1 4 0.0 4 1 4 Abs: 1 2 0.0 2 1 2 Integrating a parabolic function f(x) = alpha * x**2 Assumes that all steps have independent uncertainties! (Interval) estimate of integral: 89460.9495 89640.0505 (Plain) estimate of integral: 89550.5 89550.5