/* interpolate.c -- interpolation functions */ /* * * interpolate_local - ... (local) ** interpolate_local_order1 - ... (order 1) (linear) (2-point) ** interpolate_local_order2 - ... (order 2) (3-point) ** interpolate_local_order3 - ... (order 3) (4-point) ** interpolate_local_order4 - ... (order 4) (5-point) ** interpolate_local_order5 - ... (order 5) (6-point) */ /* To make porting easy, only the interpolate_local* functions * have been retained in this file. JT's original code (kept in * the archive directory) also has functions for derivatives and * interpolation on non-uniform grids. ---Sai. */ #include "cctk.h" #include "math.h" #define PRIVATE #define PUBLIC CCTK_REAL interpolate_local(int order, CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x); /* prototypes for private functions defined in this file */ PRIVATE CCTK_REAL interpolate_local_order1(CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x); PRIVATE CCTK_REAL interpolate_local_order2(CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x); PRIVATE CCTK_REAL interpolate_local_order3(CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x); PRIVATE CCTK_REAL interpolate_local_order4(CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x); PRIVATE CCTK_REAL interpolate_local_order5(CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x); PRIVATE CCTK_REAL interpolate_eno(CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x); /*****************************************************************************/ /* * This function does local Lagrange polynomial interpolation on a * uniform grid. That is, given order+1 uniformly spaced data points * (x0 + i*dx, y[i]) , i = 0...order, it interpolates a degree-(order) * (Lagrange) polynomial through them and evaluates this polynomial at * a specified x coordinate. In general the error in this computed * function value is O(h^(order+1)) . * * Except for possible end effects (see end_action (below)), the local * interpolation is internally centered to improve its conditioning. Even * so, the interpolation is highly ill-conditioned for small dx and/or * large order and/or x outside the domain of the data points. * * The interpolation formulas were (are) all derived via Maple, see * "./interpolate.in" and "./interpolate.out". * * Arguments: * order = (in) The order of the local interpolation, i.e. the degree * of the interpolated polynomial. * x0 = (in) The x coordinate corresponding to y[0]. * dx = (in) The x spacing between the data points. * y[0...order] = (in) The y data array. * x = (in) The x coordinate for the interpolation. */ PUBLIC CCTK_REAL interpolate_local(int order, CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x) { switch (order) { case 1: return( interpolate_local_order1(x0, dx, y, x) ); break; case 2: return( interpolate_local_order2(x0, dx, y, x) ); break; case 3: return( interpolate_local_order3(x0, dx, y, x) ); break; case 4: return( interpolate_local_order4(x0, dx, y, x) ); break; case 5: return( interpolate_local_order5(x0, dx, y, x) ); break; default: CCTK_WARN(0, "Interpolation order not supported"); return(0); } } /*****************************************************************************/ /* * This function interpolates a 1st-order (linear) Lagrange polynomial * in a 2-point [0...1] data array centered about the [0.5] grid position. */ PRIVATE CCTK_REAL interpolate_local_order1(CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x) { CCTK_REAL c0 = (+1.0/2.0)*y[0] + (+1.0/2.0)*y[1]; CCTK_REAL c1 = - y[0] + + y[1]; CCTK_REAL xc = x0 + 0.5*dx; CCTK_REAL xr = (x - xc) / dx; return( c0 + xr*c1 ); } /*****************************************************************************/ /* * This function interpolates a 2nd-order (quadratic) Lagrange polynomial * in a 3-point [0...2] data array centered about the [1] grid position. */ PRIVATE CCTK_REAL interpolate_local_order2(CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x) { CCTK_REAL c0 = y[1]; CCTK_REAL c1 = (-1.0/2.0)*y[0] + (+1.0/2.0)*y[2]; CCTK_REAL c2 = (+1.0/2.0)*y[0] - y[1] + (+1.0/2.0)*y[2]; CCTK_REAL xc = x0 + 1.0*dx; CCTK_REAL xr = (x - xc) / dx; return( c0 + xr*(c1 + xr*c2) ); } /*****************************************************************************/ /* * This function interpolates a 3rd-order (cubic) Lagrange polynomial * in a 4-point [0...3] data array centered about the [1.5] grid position. */ PRIVATE CCTK_REAL interpolate_local_order3(CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x) { CCTK_REAL c0 = (-1.0/16.0)*y[0] + (+9.0/16.0)*y[1] + (+9.0/16.0)*y[2] + (-1.0/16.0)*y[3]; CCTK_REAL c1 = (+1.0/24.0)*y[0] + (-9.0/ 8.0)*y[1] + (+9.0/ 8.0)*y[2] + (-1.0/24.0)*y[3]; CCTK_REAL c2 = (+1.0/ 4.0)*y[0] + (-1.0/ 4.0)*y[1] + (-1.0/ 4.0)*y[2] + (+1.0/ 4.0)*y[3]; CCTK_REAL c3 = (-1.0/ 6.0)*y[0] + (+1.0/ 2.0)*y[1] + (-1.0/ 2.0)*y[2] + ( 1.0/ 6.0)*y[3]; CCTK_REAL xc = x0 + 1.5*dx; CCTK_REAL xr = (x - xc) / dx; return( c0 + xr*(c1 + xr*(c2 + xr*c3)) ); } /*****************************************************************************/ /* * This function interpolates a 4th-order (quartic) Lagrange polynomial * in a 5-point [0...4] data array centered about the [2] grid position. */ PRIVATE CCTK_REAL interpolate_local_order4(CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x) { CCTK_REAL c0 = y[2]; CCTK_REAL c1 = (+1.0/12.0)*y[0] + (-2.0/ 3.0)*y[1] + (+2.0/ 3.0)*y[3] + (-1.0/12.0)*y[4]; CCTK_REAL c2 = + (-1.0/24.0)*y[0] + (+2.0/ 3.0)*y[1] + (-5.0/ 4.0)*y[2] + (+2.0/ 3.0)*y[3] + (-1.0/24.0)*y[4]; CCTK_REAL c3 = + (-1.0/12.0)*y[0] + (+1.0/ 6.0)*y[1] + (-1.0/ 6.0)*y[3] + (+1.0/12.0)*y[4]; CCTK_REAL c4 = + (+1.0/24.0)*y[0] + (-1.0/ 6.0)*y[1] + (+1.0/ 4.0)*y[2] + (-1.0/ 6.0)*y[3] + (+1.0/24.0)*y[4]; CCTK_REAL xc = x0 + 2.0*dx; CCTK_REAL xr = (x - xc) / dx; return( c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*c4))) ); } /*****************************************************************************/ /* * This function interpolates a 5th-order (quintic) Lagrange polynomial * in a 6-point [0...5] data array centered about the [2.5] grid position. */ PRIVATE CCTK_REAL interpolate_local_order5(CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x) { CCTK_REAL c0 = (+ 3.0/256.0)*y[0] + (-25.0/256.0)*y[1] + (+75.0/128.0)*y[2] + (+75.0/128.0)*y[3] + (-25.0/256.0)*y[4] + (+ 3.0/256.0)*y[5]; CCTK_REAL c1 = (- 3.0/640.0)*y[0] + (+25.0/384.0)*y[1] + (-75.0/ 64.0)*y[2] + (+75.0/ 64.0)*y[3] + (-25.0/384.0)*y[4] + (+ 3.0/640.0)*y[5]; CCTK_REAL c2 = (- 5.0/ 96.0)*y[0] + (+13.0/ 32.0)*y[1] + (-17.0/ 48.0)*y[2] + (-17.0/ 48.0)*y[3] + (+13.0/ 32.0)*y[4] + (- 5.0/ 96.0)*y[5]; CCTK_REAL c3 = (+ 1.0/ 48.0)*y[0] + (-13.0/ 48.0)*y[1] + (+17.0/ 24.0)*y[2] + (-17.0/ 24.0)*y[3] + (+13.0/ 48.0)*y[4] + (- 1.0/ 48.0)*y[5]; CCTK_REAL c4 = (+ 1.0/ 48.0)*y[0] + (- 1.0/ 16.0)*y[1] + (+ 1.0/ 24.0)*y[2] + (+ 1.0/ 24.0)*y[3] + (- 1.0/ 16.0)*y[4] + (+ 1.0/ 48.0)*y[5]; CCTK_REAL c5 = (- 1.0/120.0)*y[0] + (+ 1.0/ 24.0)*y[1] + (- 1.0/ 12.0)*y[2] + (+ 1.0/ 12.0)*y[3] + (- 1.0/ 24.0)*y[4] + (+ 1.0/120.0)*y[5]; CCTK_REAL xc = x0 + 2.5*dx; CCTK_REAL xr = (x - xc) / dx; return( c0 + xr*(c1 + xr*(c2 + xr*(c3 + xr*(c4 + xr*c5)))) ); } /*****************************************************************************/ CCTK_REAL interpolate_eno(CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x) { CCTK_REAL diff[4]; CCTK_INT seed = 0; CCTK_INT j,r; CCTK_REAL c0; CCTK_REAL c1; CCTK_REAL c2; CCTK_REAL xc; CCTK_REAL xr; /* Find seed index */ while (x > x0 + dx * ((float)seed+0.5)) seed++; if (seed!=2) { /* Not enough stencil, only perform linear interpolation */ seed = 0; while (x > x0 + dx * ((float)(seed+1))) seed++; if (seed==4) seed=3; return y[seed] + (y[seed+1]-y[seed])/dx * (x-x0-(float)seed*dx); } /* Calculate the undivided differences */ for (j=0; j<=3; j++) diff[j] = y[j+1] - y[j]; /* Find the stencil */ if ( fabs(diff[1]) < fabs(diff[2]) ) { if ( fabs(diff[1]-diff[0]) < fabs(diff[2]-diff[1]) ) r = 0; else r = 1; } else { if ( fabs(diff[2]-diff[1]) < fabs(diff[3]-diff[2]) ) r = 1; else r = 2; } /* Interpolate second order */ c0 = y[r+1]; c1 = (-1.0/2.0)*y[r] + (+1.0/2.0)*y[r+2]; c2 = (+1.0/2.0)*y[r] - y[r+1] + (+1.0/2.0)*y[r+2]; xc = x0 + dx * (CCTK_REAL)(r+1); xr = (x - xc) / dx; return( c0 + xr*(c1 + xr*c2) ); }