Interpolación cúbica monótona

En el campo matemático del análisis numérico, la interpolación cúbica monótona es una variante del interpolador cúbico de Hermite que preserva la monotonicidad del conjunto de datos que se está interpolando.​

Ejemplo que muestra la interpolación cúbica no monótona (en rojo) y la interpolación cúbica monótona (en azul) de un conjunto de datos monótonos

La interpolación lineal preserva la monotonicidad (es decir, que la altura de la serie de puntos y de la curva sea siempre ascendente o descendente, sin cambios en el signo de la pendiente), pero la interpolación cúbica no la garantiza.

Interpolación cúbica monótona de Hermite

La interpolación monótona se puede lograr usando el interpolador cúbico de Hermite con las tangentes $m_{i}$ modificadas para garantizar la monotonicidad del spline de Hermite resultante.

También está disponible un algoritmo para la interpolación monótona quíntica de Hermite.

Selección del polinomio interpolante

Hay varias formas de seleccionar tangentes de interpolación para cada punto del conjunto de datos. Esta sección describirá el uso del método de Fritsch-Carlson.​ Téngase en cuenta que solo se requiere una pasada del algoritmo.

Considérese que los puntos del conjunto de datos $(x_{k},y_{k})$ se indexen en orden para $k=1,\,\dots \,n$.

1. Calcular las pendientes de las rectas secantes entre puntos sucesivos:

   $\delta _{k}={\frac {y_{k+1}-y_{k}}{x_{k+1}-x_{k}}}$

   para $k=1,\,\dots \,n-1$.
   
   
2. Estas asignaciones son provisionales, y pueden ser reemplazadas en los pasos restantes. Inicialícense las tangentes en cada punto de datos interior como el promedio de las secantes,

   $m_{k}={\frac {\delta _{k-1}+\delta _{k}}{2}}$

   para $k=2,\,\dots \,n-1$.
   
   Para los puntos finales, utilizar las diferencias de un solo lado:

   $m_{1}=\delta _{1}\quad {\text{ y }}\quad m_{n}=\delta _{n-1}\,$.

   Si $\delta _{k-1}$ y $\delta _{k}$ tienen signos opuestos, configurar $m_{k}=0$.
   
   
3. Para $k=1,\,\dots \,n-1$, siempre que $\delta _{k}=0$ (donde dos $y_{k}=y_{k+1}$ sucesivos sean iguales),
   establecer $m_{k}=m_{k+1}=0,$ ya que el spline que conecta estos puntos siempre debe decrecer (o crecer) a medida que se avanza según el eje x para preservar la monotonicidad.
   Ignorar los pasos 4 y 5 para aquellos $k\,$.
   
   
4. Dejar

   $\alpha _{k}=m_{k}/\delta _{k}\quad {\text{ y }}\quad \beta _{k}=m_{k+1}/\delta _{k}$.

   Si $\alpha _{k}$ o $\beta _{k}$ es negativo, entonces los puntos de datos de entrada no son estrictamente monótonos y $(x_{k},\,y_{k})$ es un extremo local. En tales casos, aún se pueden generar curvas monótonas "por partes" eligiendo $m_{k}=0\,$ si $\alpha _{k}<0$ o $m_{k+1}=0\,$ si $\beta _{k}<0$, aunque la monotonicidad estricta no es posible globalmente.
   
   
5. Para evitar sobrepaso y garantizar la monotonicidad, se debe cumplir al menos una de las tres condiciones siguientes:

(a) La función

$\phi _{k}=\alpha _{k}-{\frac {(2\alpha _{k}+\beta _{k}-3)^{2}}{3(\alpha _{k}+\beta _{k}-2)}}>0\,$, 'o

(b) $\alpha _{k}+2\beta _{k}-3\leq 0\,$, 'o

(c) $2\alpha _{k}+\beta _{k}-3\leq 0\,$.

Solo la condición (a) es suficiente para garantizar una monotonicidad estricta: $\phi _{k}$ debe ser positivo.

Una forma sencilla de satisfacer esta restricción es restringir el vector $(\alpha _{k},\,\beta _{k})$ a un círculo de radio 3. Es decir, si es $\alpha _{k}^{2}+\beta _{k}^{2}>9\,$, establecer

$\tau _{k}={\frac {3}{\sqrt {\alpha _{k}^{2}+\beta _{k}^{2}}}}\,$,

y cambiar la escala de las tangentes mediante

$m_{k}=\tau _{k}\,\alpha _{k}\,\delta _{k}\quad {\text{ y }}\quad m_{k+1}=\tau _{k}\,\beta _{k}\,\delta _{k}\,$.

Como alternativa, es suficiente restringir $\alpha _{k}\leq 3$ y $\beta _{k}\leq 3\,$. Para lograr esto, si es $\alpha _{k}>3{\text{ o }}\beta _{k}>3\,$, configurar $m_{k}=3\,\delta _{k}\,$.

Interpolación cúbica

Después del preprocesamiento anterior, la evaluación del spline interpolado es equivalente a la del interpolador cúbico de Hermite, utilizando los datos $x_{k}$, $y_{k}$ y $m_{k}$ para $k=1,\,\dots \,n$.

Para evaluar en $x$, buscar el índice $k$ en la secuencia donde $x$ se encuentra entre $x_{k}$ y $x_{k+1}$, es decir: $x_{k}\leq x\leq x_{k+1}$. Calcular

$\Delta =x_{k+1}-x_{k}\quad {\text{ y }}\quad t={\frac {x-x_{k}}{\Delta }}$

entonces, el valor interpolado es

$f_{\text{ interpolada }}(x)=y_{k}\cdot h_{00}(t)+\Delta \cdot m_{k}\cdot h_{10}(t)+y_{k+1}\cdot h_{01}(t)+\Delta \cdot m_{k+1}\cdot h_{11}(t)$

donde $h_{ii}$ son las funciones básicas para el interpolador cúbico de Hermite.

Código de ejemplo

El siguiente código en JavaScript toma un conjunto de datos y produce una función de interpolación monótona mediante un spline cúbico:

/*
 * Monotone cubic spline interpolation
 * Usage example listed at bottom; this is a fully-functional package. For
 * example, this can be executed either at sites like
 * https://www.programiz.com/javascript/online-compiler/
 * or using nodeJS.
 */
function DEBUG(s) {
    /* Uncomment the following to enable verbose output of the solver: */
    //console.log(s);
}
var outputCounter= 0;
var createInterpolant= function(xs, ys) {
    var i, length= xs.length;
    
    // Deal with length issues
    if (length != ys.length) {throw 'Need an equal count of xs and ys.'; }
    if (length=== 0) {return function(x) {return 0; }; }
    if (length=== 1) {
        // Impl: Precomputing the result prevents problems if ys is mutated later and allows garbage collection of ys
        // Impl: Unary plus properly converts values to numbers
        var result= +ys[0];
        return function(x) {return result; };
    }
    
    // Rearrange xs and ys so that xs is sorted
    var indexes= [];
    for (i= 0; i < length; i++) {indexes.push(i); }
    indexes.sort(function(a, b) {return xs[a] < xs[b] ? -1 : 1; });
    var oldXs= xs, oldYs= ys;
    // Impl: Creating new arrays also prevents problems if the input arrays are mutated later
    xs= []; ys= [];
    // Impl: Unary plus properly converts values to numbers
    for (i= 0; i < length; i++) {
        xs[i]= +oldXs[indexes[i]];
        ys[i]= +oldYs[indexes[i]];
    }

    DEBUG("debug: xs= [ " + xs + " ]")
    DEBUG("debug: ys= [ " + ys + " ]")
    
    // Get consecutive differences and slopes
    var dys= [], dxs= [], ms= [];
    for (i= 0; i < length - 1; i++) {
        var dx= xs[i + 1] - xs[i], dy= ys[i + 1] - ys[i];
        dxs[i]= dx;
        dys[i]= dy;
        ms[i]= dy/dx;
    }
    // Get degree-1 coefficients
    var c1s= [ms[0]];
    for (i= 0; i < dxs.length - 1; i++) {
        var m= ms[i], mNext= ms[i + 1];
        if (m*mNext <= 0) {
            c1s[i]= 0;
        } else {
            var dx_= dxs[i], dxNext= dxs[i + 1], common= dx_ + dxNext;
            c1s[i]= 3*common/((common + dxNext)/m + (common + dx_)/mNext);
        }
    }
    c1s.push(ms[ms.length - 1]);

    DEBUG("debug: dxs= [ " + dxs + " ]")
    DEBUG("debug: ms= [ " + ms + " ]")
    DEBUG("debug: c1s.length= " + c1s.length)
    DEBUG("debug: c1s= [ " + c1s + " ]")
    
    // Get degree-2 and degree-3 coefficients
    var c2s= [], c3s= [];
    for (i= 0; i < c1s.length - 1; i++) {
        var c1= c1s[i];
        var m_= ms[i];
        var invDx= 1/dxs[i];
        var common_= c1 + c1s[i + 1] - m_ - m_;
        DEBUG("debug: " + i + ". c1= " + c1);
        DEBUG("debug: " + i + ". m_= " + m_);
        DEBUG("debug: " + i + ". invDx= " + invDx);
        DEBUG("debug: " + i + ". common_= " + common_);
        c2s[i]= (m_ - c1 - common_)*invDx;
        c3s[i]= common_*invDx*invDx;
    }
    DEBUG("debug: c2s= [ " + c2s + " ]")
    DEBUG("debug: c3s= [ " + c3s + " ]")

    // Return interpolant function
    return function(x) {
        // The rightmost point in the dataset should give an exact result
        var i= xs.length - 1;
        // Uncomment the following to return only the interpolated value.
        //if (x== xs[i]) {return ys[i]; }
        
        // Search for the interval x is in, returning the corresponding y if x is one of the original xs
        var low= 0, mid, high= c3s.length - 1;
        while (low <= high) {
            mid= Math.floor(0.5*(low + high));
            var xHere= xs[mid];
            if (xHere < x) {low= mid + 1; }
            else if (xHere > x) {high= mid - 1; }
            else {
                // Uncomment the following to return only the interpolated value.
                //return ys[mid];
                low= c3s.length - 1;
                high= mid;
                break;
            }
        }
        i= Math.max(0, high);

        // Interpolate
        var diff= x - xs[i];
        outputCounter++;
        var interpolatedValue= ys[i] + diff * (c1s[i] + diff *  (c2s[i] + diff * c3s[i]));
        // The value of the interpolator's derivative at this point.
        var derivativeValue= c1s[i] + diff * (2*c2s[i] + diff * 3*c3s[i]);
        DEBUG("debug: #" + outputCounter + ". x= " + x + ". i= " + i + ", diff= " + diff + ", interpolatedValue= " + interpolatedValue + ", derivativeValue= " + derivativeValue);
        // Uncomment the following to return only the interpolated value.
        // return interpolatedValue;
        return [ interpolatedValue, derivativeValue ];
    };
};

/*
   Usage example below will approximate x^2 for 0 <= x <= 4.

   Command line usage example (requires installation of nodejs):
   node monotone-cubic-spline.js
*/

var X= [0, 1, 2, 3, 4];
var F= [0, 1, 4, 9, 16];
var f= createInterpolant(X,F);
var N= X.length;
console.log("# BLOCK 0 :: Data for monotone-cubic-spline.js");
console.log("X" + "\t" + "F");
for (var i= 0; i < N; i += 1) {
    console.log(X[i] + '\t' + F[i]);
}
console.log(" ");
console.log(" ");
console.log("# BLOCK 1 :: Interpolated data for monotone-cubic-spline.js");
console.log("      x       " + "\t\t" + "     P(x)      " + "\t\t" + "    dP(x)/dx     ");
var message= '';
var M= 25;
for (var i= 0; i <= M; i += 1) {
    var x= X[0] + (X[N-1]-X[0])*i/M;
    var rvals= f(x);
    var P= rvals[0];
    var D= rvals[1];
    message += x.toPrecision(15) + '\t' + P.toPrecision(15) + '\t' + D.toPrecision(15) + '\n';
}
console.log(message);

Bibliografía

• Fritsch, F. N.; Carlson, R. E. (1980). «Monotone Piecewise Cubic Interpolation». SIAM Journal on Numerical Analysis (SIAM) 17 (2): 238-246. doi:10.1137/0717021. 
• Dougherty, R.L.; Edelman, A.; Hyman, J.M. (April 1989). «Positivity-, monotonicity-, or convexity-preserving cubic and quintic Hermite interpolation». Mathematics of Computation 52 (186): 471-494. doi:10.2307/2008477. 

Enlaces externos

• Implementación GNU General Public License con licencia C++v2: MonotCubicInterpolator.cpp -origins/blob/master/dune-cornerpoint/common/MonotCubicInterpolator.hpp MonotCubicInterpolator.hpp

• Datos: Q4044484