/* interpolate.c -- interpolation functions */ /* $Header$ */ /* * 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 #include #include "cctk.h" CCTK_REAL interpolate_local(int order, CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x); /*****************************************************************************/ /* * This function interpolates a 1st-order (linear) Lagrange polynomial * in a 2-point [0...1] data array centered about the [0.5] grid position. */ static 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. */ static 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. */ static 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. */ static 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. */ static 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)))); } /* * 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. */ CCTK_REAL interpolate_local(int order, CCTK_REAL x0, CCTK_REAL dx, CCTK_REAL y[], CCTK_REAL x) { int start, end; start = 0; end = order; /* The possible interpolation order that can be used */ if (x0 > x) return y[start]; if (dx * end + x0 < x) return y[end]; switch (order) { case 1: return interpolate_local_order1(x0, dx, y, x); case 2: return interpolate_local_order2(x0, dx, y, x); case 3: return interpolate_local_order3(x0, dx, y, x); case 4: return interpolate_local_order4(x0, dx, y, x); case 5: return interpolate_local_order5(x0, dx, y, x); default: CCTK_WARN(0, "Interpolation order not supported"); return 0; /* Avoid warning */ } }