Anonymous user
Thiele's interpolation formula: Difference between revisions
→{{header|D}}: added D
m (→{{header|Tcl}}: closer, maybe…) |
(→{{header|D}}: added D) |
||
Line 77:
pi estimate using tan interpolation: 3.1415926535898
</pre>
=={{header|D}}==
<lang d>import std.stdio ;
import std.math ;
U[] myMap(U,V) (U delegate(V) f, V[] a) {
V[] r = new V[](a.length) ;
foreach(i, v ; a)
r[i] = f(v) ;
return r ;
}
T[] myRng(T)(T start, T end, T step) {
T[] r ;
for(; start < end ; start += step)
r ~= start ;
return r ;
}
alias real delegate(real) RealFun ;
const real End = 1.55L , Step = 0.05L , Start = 0.0L ;
class Thiele {
const real[] F, X, RhoY, RhoX ;
int NN ;
this(real[] f, real[] x) {
F = f.dup ;
X = x.dup ;
NN = X.length ;
assert(NN > 3, "at leat 4 values") ;
assert(X.length == F.length, "input arrays not of same size") ;
RhoY = rhoN(F, X) ;
RhoX = rhoN(X, F) ;
}
this(RealFun fun, real start = Start, real end = End, real step = Step) {
this(myMap(fun, myRng(start, end, step)),
myRng(start, end, step)) ;
}
private static real[] rhoN(real[] f, real[] x) {
int N = x.length;
real[][] p = new real[][] (N, N) ;
for(int i = 0; i < N ; i++)
p[i][0] = f[i] ;
for(int i = 0; i < N - 1 ; i++)
p[i][1] = (x[i] - x[i+1]) / (p[i][0] - p[i+1][0]) ;
for(int j = 2; j < N - 1; j++)
for(int i = 0; i < N - j - 1 ; i++)
p[i][j] = p[i+1][j-2] + (x[i] - x[i+j]) /
(p[i][j-1] - p[i+1][j-1]) ; ;
return p[1].dup ;
}
real evalY(real x) {
real a = 0.0L ;
for(int i = NN - 3; i >= 2 ; i--)
a = (x - X[i]) / (RhoY[i] - RhoY[i-2] + a) ;
return F[1] + (x - X[1]) / (RhoY[1] + a) ;
}
real evalX(real y) { // inverse
real a = 0.0L ;
for(int i = NN - 3; i >= 2 ; i--)
a = (y - F[i]) / (RhoX[i] - RhoX[i-2] + a) ;
return X[1] + (y - F[1]) / (RhoX[1] + a) ;
}
}
void main() {
auto fun = (real x) { return std.math.sin(x) ; } ;
auto t = new Thiele(fun) ;
writefln(" %d interpolating points", t.NN) ;
writefln("std.math.sin(0.5) : %20.18f", std.math.sin(0.5L)) ;
writefln(" Thiele sin(0.5) : %20.18f", t.evalY(0.5L)) ;
writefln("*%20.18f library constant", std.math.PI) ;
writefln(" %20.18f 6*sin`(0.5)", t.evalX(0.5L) * 6.0L) ;
t = new Thiele((real x) {return std.math.cos(x) ; }) ;
writefln(" %20.18f 3*cos`(0.5)", t.evalX(0.5L) * 3.0L) ;
t = new Thiele((real x) {return std.math.tan(x) ; }) ;
writefln(" %20.18f 4*tan`(0.5)", t.evalX(1.0L) * 4.0L) ;
}</lang>
output:
<pre> 32 interpolating points
std.math.sin(0.5) : 0.479425538604203000
Thiele sin(0.5) : 0.479425538604203000
*3.141592653589793239 library constant
3.141592653589793238 6*sin`(0.5)
3.141592653589793238 3*cos`(0.5)
3.141592653589793238 4*tan`(0.5)</pre>
=={{header|Tcl}}==
| |||