Rosetta Code

Thiele's interpolation formula: Difference between revisions

Page Discussion

Thiele's interpolation formula (view source)

Revision as of 14:46, 3 October 2010

225 bytes removed ,  13 years ago
→‎{{header|Tcl}}: Corrected version
(Reference of reciprocal difference.)
(→‎{{header|Tcl}}: Corrected version)
Line 246:
 
=={{header|Tcl}}==
{{works with|Tcl|8.5}}<br>{{trans|D}}▼
{{incorrect|Tcl|It produces the wrong values for output.}}
<lang tcl>#
▲{{works with|Tcl|8.5}}
### Create a thiele-interpretation function with the given name that interpolates
<lang tcl># Multi-list map; deeper voodoo than normal!
### off the given table.
proc map args {
#
foreach {a b} [lrange $args 0 end-1] {▼
proc thiele {name : X Y-> xF} {▼
lappend parms [incr v] $b
upvar 1 $a $v
}▼
foreach {*}$parms {lappend result [uplevel 1 [lindex $args end]]}
return $result
}▼
 
▲proc thiele {X Y x} {
# Sanity check
if {[llength $X] != [llength $YF]} {
error "unequal length lists supplied: [llength $X] != [llength $YF]"
}
 
Line 267 ⟶ 260:
### Compute the table of reciprocal differences
#
set nN [expr {[llength $Y] - 1}X]
set p [lrepeat $N [lrepeat $N 0.0]]
# Zero-th order
for {set ri [list0} {$Y]i<$N} {incr i} {
lset p $i 0 [lindex $F $i]
# First order
▲ }
lappend r [map a [lrange $X 0 end-1] b [lrange $X 1 end] \
c [lrangefor $Y{set i 0} end{$i<$N-1]} d [lrange $Y 1{incr end]i} \{
lset p $i e1 [lrepeat $n 0]expr {
expr {($a - $b) /([lindex ($cX $i] - [lindex $dX $i+1]) $e)}/
set a [expr {($x - $b) / ([lindex $cp $i 0] - [lindex $dp $i+1 0] + $a)}]▼
}]▼
# Higher order
for {set i 2} {$i <= $n} {incr i} {▼
set p [lindex $r end]; # previous row
set pp [lindex $r end-1]; # doubly-previous row; longer than $p by 1
lappend r [map a [lrange $X 0 end-$i] b [lrange $X $i end] \
c [lrange $p 0 end-1] d [lrange $p 1 end] \
e [lrange $pp 1 end-1] {▼
expr {($a - $b) / ($c - $d + $e)}
}]
}▼
▲ for {set ij 2} {$i j<= $nN-1} {incr ij} {
for {set i 0} {$i<$N-$j} {incr i} {
lset p $i $j [expr {
[lindex $p $i+1 $j-2] +
([lindex $X $i] - [lindex $X $i+$j]) /
([lindex $p $i $j-1] - [lindex $p $i+1 $j-1])
▲ }]
▲ }
}
 
#
### ActuallyMake evaluatepseudo-curried function that actually evaluates Thiele's formula
#
interp alias {} $name {} apply {{X rho f1 x} {
set a 0.0
foreach b [lreverse [lrange $X 2 end]] c [lreverse [lrange $r 2 end]] \
foreach Xi d [lreverse [lrange $rX 02 end-2]] {\
▲ Ri e[lreverse [lrange $pprho 12 end-1]] {\
▲ set a [expr {($x - $b) / ([lindex $c 0] - [lindex $d 0] + $a)}]
▲ Ri2 foreach {a b}[lreverse [lrange $argsrho 0 end-12]] {
▲ }
set a [expr {[lindex $Y 0] + ($x - [lindex $X 0]Xi) / ([lindex $rRi 1- 0]$Ri2 + $a)}]
}
expr {$f1 + ($x - [lindex $X 1]) / ([lindex $rho 1] + $a)}
}} $X [lindex $p 1] [lindex $F 1]
}</lang>
Demonstration code:
Line 306 ⟶ 303:
}
 
interp alias {}thiele invSin {} thiele: $trigTable(sin) -> $trigTable(x)
interp alias {}thiele invCos {} thiele: $trigTable(cos) -> $trigTable(x)
interp alias {}thiele invTan {} thiele: $trigTable(tan) -> $trigTable(x)
}
initThieleTest
Line 316 ⟶ 313:
Output:
<pre>
pi estimate using sin interpolation: 3.76073733417127981415926535897936
pi estimate using cos interpolation: 13.3300617930864258141592653589793
pi estimate using tan interpolation: 83.98541899822322141592653589794
</pre>
Anonymous user
Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/158164"