From: spellucci@mathematik.tu-darmstadt.de (Peter Spellucci)
Newsgroups: sci.math.num-analysis,sci.math
Subject: Re: Code for N-Dimensional linear interpolation
Date: 9 Dec 1998 17:53:03 GMT

In article <74ks5i$m45$1@news-1.news.gte.net>,
 Ken Van de Houten <kjvand@gte.net> writes:
|> 
|> Does anyone have a C or FORTRAN routine for performing
|> a linear interpolation over N dimensional data sampled at unequal
|> spacings ?      i.e. I have some N-Dimensional arrays of data, with
|> different limits and sample rates along each axis.    I need to compute
|> an interpolated value within the space.
|> 
|> Linear interpolation is certainly the simplest,
no, this is not true in your case.
what you obviously have are data of the form (x_i,y_j,z_k,...) with
x_i form an interval [a1,b1], y_j from an interval [a2,b2], and so on,
that means data from a grid on a hypercube. linear interpolation requires N+1
data points from a simplex in R^N, for N=2 a triangle, N=3 a tetrahedron and so
on. hence you would first have to generate a simplicial grid from your
hypercube-grid. this is of course possible but not completely trivial. 
much easier is the socalled multilinear interpolation. here the 
function is a product of a linear function in x times a linear function in y
times ... 
let x_i <= x <= x_{i+1}, y_j <= y <= y_{j+1} , z_k <= z <= z_{k+1}
i.e. the case N=3. then the interpolation function becomes 
    
  f(x_i,y_j,z_k)*(x-x_{i+1})*(y-y_{j+1})*(z-z_{k+1})/
               ((x_i-x_{i+1})*(y_j-y_{j+1})*(z_k-z_{k+1}))
 +
  f(x_{i+1},y_j,z_k)*(x-x_i)*(y-y_{j+1})*(z-z_{k+1})/
               ((x_{i+1}-x_i)*(y_j-y_{j+1})*(z_k-z_{k+1}))
 +
  f(x_i,y_{j+1},z_k)*(x-x_{i+1})*(y-y_j)*(z-z_{k+1})/
               ((x_i-x_{i+1})*(y_{j+1}-y_j)*(z_k-z_{k+1}))
 +

  f(x_{i+1},y_{j+1},z_k)*(x-x_i)*(y-y_j)*(z-z_{k+1})/
               ((x_i-x_{i+1})*(y_j-y_{j+1})*(z_k-z_{k+1})
 +

  f(x_i,y_j,z_{k+1})*(x-x_{i+1})*(y-y_{j+1})*(z-z_k)/
               ((x_i-x_{i+1})*(y_j-y_{j+1})*(z_{k+1}-z_k))
 +
  f(x_{i+1},y_j,z_{k+1})*(x-x_i)*(y-y_{j+1})*(z-z_k)/
               ((x_{i+1}-x_i)*(y_j-y_{j+1})*(z_{k+1}-z_k))
 +
 +
  f(x_i,y_{j+1},z_{k+1})*(x-x_{i+1})*(y-y_j)*(z-z_k)/
               ((x_i-x_{i+1})*(y_{j+1}-y_j)*(z_{k+1}-z_k))
 +

  f(x_{i+1},y_{j+1},z_{k+1})*(x-x_i)*(y-y_j)*(z-z_k)/
               ((x_i-x_{i+1})*(y_j-y_{j+1})*(z_{k+1}-z_k)
the structure of that formula is completely analogous to the Lagrange
interpolation formula in one variable: each of the 8 (in general 2^N))
functions values is multiplied by a polynomial linear in each variable which is 
1 at the corresponding grid point and zero at the other grid points on the
hypercube. the interpolation formula is valid on each hypercube, but 
changes with i,j,k of course. 
this can be extended to higher order interpolation in a similar way.
hope this helps
peter