Interpolación

Contexto Suponemos que tenemos una función dada por una tabla, donde se indican un conjunto de valores (xi,yi). El problema consiste en estimar qué valores puede tomar Y para x que no coinciden con los de la tabla. Si xi<x<xi+1, esta estimación es a lo que llamamos interpolación. Si está fuera del intervalo, se llama extrapolación.
Interpolación con polinomios de Lagrange Consiste en determinar el polinomio que pasa por los puntos de la tabla. Este polinomio se puede determinar planteando un sistema de ecuaciones, o bien por una técnica de Lagrange. Ver aquí.
Librería En interpolacion.js está desarrollado el código para calcular el polinomio de Lagrange.
Se utiliza la librería ../Polinomios/polinomios.js que debe incluirse en el HTML.
Como datos de entrada debemos introducir el array con los datos, de la siguiente manera: 
[[x0,y0],[x1,y1],[x2,y2],..,[x(n-1), y(n-1)]]

Se utiliza la función 
function interpolacionLagrange(datos)
que devuelve un array unidimensional con los coeficientes del polinomio resultante, que podemos escribir con la función polinomioToString(polinomio) y evaluarlo con evalPol(polinomio,x).
Ver siguiente ejemplo: 

let salida = "";
let datos = [[-4,10],[-3,0],[-2,0],[-1,0],[1,0],[2,0],[3,0],[4,-10]];
for (let dato of datos) salida += "P(" + dato[0] + ") = " + dato[1] + "\n"
let polLagrange = interpolacionLagrange(datos);
salida += "\n" + polinomioToString(polLagrange) + "\n";
for (let x = -5; x < 5; x++) salida += "P(" + x + ") = " + evalPol(polLagrange, x) + "\n";
console.log(salida);
Interpolación con splines cúbicos
Desarrollo Se puede ver un tratado de splines aquí.
Tenemos n+1 puntos (de 0 a n): 
[[x0,y0],[x1,y1],[x2,y2],..,[x(n), y(n)]]

Queremos determinar una función a trozos, de forma que para cada intervalo xi < x < xi+1 tengamos por un polinomio cúbico de interpolación para los x de ese intervalo. Para n+1 puntos tendremos n intervalos, y por tanto, polinomios distintos. El polinomio del intervalo i lo denotamos: \[Y_i(x) = A_i·x^3 + B_i·x^2 + C_i·x + D_i\]
Tenemos 4·n incógnitas.
Podemos eliminar una incógnita de cada polinomio escribiendo: \[Y_i(x) = a_i·(x-x_i)^3 + b_i·(x-x_i)^2 + c_i·(x-x_i) + d_i\]
El primer polinomio se puede obtener del segundo mediante las siguientes expresiones:
Ai = ai
Bi = bi-3·ai·xi
Bi = 3·xi2·ai - 2·bi·xi+ci
Di = -ai·xi3 + bi·xi2-c·xi+di
Poniendo como condición que el polinomio contenga el punto (xi, yi) obtenemos
Yi(xi) = yi = di, ∀i ϵ [0,..,n-1] .
Llamamos hi = xi+1-xi
La función y sus derivadas en el intervalo son:
Yi(x) = ai·(x-xi)3 + bi·(x-xi)2 + ci·x + di
Y'i(x) = 3·ai·(x-xi)2 + 2·bi·(x-xi) + ci
Y"i(x) = 6·ai·(x-xi) + 2·bi.
Ponemos la condición de que Yi(xi+1) = yi+1 ⇒ yi+1 = ai·hi3+bihi2+c·hi+di
Llamamos si = Y"i(xi) = 2·bi
Ponemos la condición de que la segunda derivada sea continua: si+1=6·ai·hi+2·bi.
De esta dos condiciones deducimos:
bi = si/2
ai=(si+1-si)/(6·hi).
Sustituyendo estas expresiones, obtenemos ci=(yi+1-yi)/hi-(2·hi·si+hi·si+1)/6
Ponemos la última condición: Y'i-1(xi) = Y'i(xi)
Operando se obtiene el siguiente sistema de ecuaciones:
hi-1·si-1 + 2·(hi-1+hi)·si + hi·si+1 = = 6·((yi+1-yi)/hi-(yi-yi-1)/hi-1 

					[A] · [s] = [b]
	
[[ h0   2·(h0+h1)      h1                                                    ]  [[  s0  ]
 [          h1     2·(h1+h2)     h2                                          ]   [  s1  ]
 [                           ...                                            ]   [  ..  ]
 [                              ...                                         ] x [  ..  ]  =
 [                                 ...                                      ]   [  ..  ]
 [                                            hn-2     2·(hn-2+hn-1)     hn-1 ]]  [  sn  ]]

 [[     6·((y2-y1)/h1 - (y1-y0)/h0)      ]
  [     6·((y3-y2)/h2 - (y2-y1)/h1)      ]
  [                 ...                 ]
  [                 ...                 ]
  [                 ...                 ]
  [ 6·((yn-yn-1)/hn-1 - (yn-1-yn-2)/hn-2) ]]
sn podría parecer que se refiere a un polinomio tras el punto (xi,yi), por lo que conviene aclarar que lo sefinimos como el valor de y"n-1(xn).
Ahora hemos obtenido un sistema de n+1 incógnitas, es decir, lo hemos reducido a la cuarta parte.
Si obtenemos los valores de si calcularemos:
ai = (si+1 - si)/(6·hi)
bi = si/2
ci = (yi+1 - yi) / hi - (2·hi·si + hi·si+1)/6
di = yi

El sistema encontrado:

  • Es tridiagonal.
  • Es simétrica con diagonal principal dominante.
  • Tiene n-1 ecuaciones y n+1 incógnitas.
Al tener menos ecuaciones que incógnitas, implica que podemos poner dos condiciones más, una por cada uno de los polinomios extremos, que pondremos en forma de ecuaciones. Y tenemos varias opciones.
  • Spline natural Esto es asignar s0 = 0; sn=0. Con ello logramos la función más suave posible, según demuestran en el documento indicado más arriba. Podemos aplicarlo a los dos extremos o aplicarlo a sólo uno de los extremos. Añadir estas ecuaciones es una fila con sólo un 1 en el coeficiente correspondiente, resto ceros, y término independiente 0.
  • Si tenemos la mitad de los datos, porque la función es simétrica y derivable, el extremo donde tenemos la simetría la derivada debe ser cero (pendiente horizontal). Si es en x0. La función derivada es
    Y'i = 3·ai·(x-xi)2+2·bi·(x-xi)+ci
    para la última función (yn-1(xn) ) sustituyendo los coeficientes en función de si y operando, tenemos que
    • Simetría por la izquierda:
      y'0(x0) = 0 = c0. Operando y'i(xi) = (yi+1-yi)/hi - hi/3·si - hi/6·si+1, particularizando para i=0 ⇒ 1/3·h0·s0 + 1/6·h0·s1 = (y1-y0)/h0
    • Simetría por la derecha:
      y'i(xi+1) = hi/3·si+1 + hi/6·si + (yi+1-yi)/hi Particurarizando para i = n-1 obtenemos:
      y'n-1(x=xn) = 0 = hn-1/3·sn + hn-1/6·sn-1 + (yn-yn-1)/hn-1
      (Sin comprobar).
  • Si la función es periódica, tomaremos el punto final igual que el inicial, y0(x0) = yn(xn) Además haremos coincidir las pendientes (y'0(x0) = y'n-1(xn)), e igualaremos las segundas derivadas (s0 = sn).
    En este caso el sistema de ecuaciones sería tridiagonal, salvo por dos ecuaciones:
    h0/3·s0 + h0/6·s1 + hn-1/6·sn-1 + hn-1/3·sn = (y1-y0)/h0 - (yn-yn-1)/hn-1
    s0 - sn = 0
    (sin comprobar)
  • Otras condiciones que se suelen usar consiste en mantener las derivadas segundas de las funciones extremas igual que las de sus contiguas: s0 = s1; sn-1 = sn-2.

Librería Tenemos como datos de entrada: 
[[x0,y0],[x1,y1],[x2,y2],..,[x(n-1), y(n-1)]]

Se ha definido la clase SplineCubico cuyo constructor es: 

let spline = new SplineCubico(datos, condicionIzda, condicionDcha);

donde:
  • datos es el array de datos de entrada.
  • condicionIzda puede ser "natural" para poner la condición y"(x0) = 0; "simetrica" que pone la condición y'(x0) = 0; "parabolica" para poner condición y"(x0] = y"(x1>)]; y "periódica" para condición y"(x0] = y"(xn-1>)] y y'(x0] = y'(xn-1>)].
  • condicionDcha puede ser "natural" para poner la condición y"(xn-1) = 0; "simetrica" que pone la condición y'(xn-1) = 0; "parabolica" para poner condición y"(xn-1] = y"(xn-2>)]; y "periódica" para condición y"(x0] = y"(xn-1>)] y y'(x0] = y'(xn-1>)].

En el caso de que sea periódica, tenemos que introducir los datos de forma que y(x0) = y(xn-1). Además, si se pone la condición "periodica" en cualquera de los extremos, se aplica también al otro, ignorándose cualquier otra condición que pudiéramos tener establecida.

Al instanciar el objeto spline se plantea y resuelve el sistema de ecuaciones, para lo que se debe incluir el fichero sistemaEcuacionesLineales.js. Si se aplica la condición de periodicidad el sistema se resuelve por Gauss Seidel en el caso general, mientras para el resto de condiciones se resuelve optimizando que la matriz del sistema de ecuaciones es tridiagonal.

Tras el cálculo del sistema de ecuaciones, el objeto spline guarda lo arrays x[], a[], b[], c[] y d[], de forma que para un x busca el intervalo i que cumple x[i] < x <= x[i+1] y devuelve a[i]*x^3+b[i]^2+c[i]*x+d[i]. Esto es lo que hace el método spline.evalua(x).

El método spline.derivada(x) devuelve la derivada en x del polinomio del splite correspondiente.
Los métodos spline.derivada2(x) y spline.derivada3(x) devuelven la segunda y tercera derivada.

El método spline.integral(xInf, xSup) devuelve el área bajo la curva entre xInf y xSup.

El método spline.longitud(xInf, xSup) devuelve la longitud de la curva entre xInf y xSup.

El método spline.radioCurv(x) devuelve el radio de curvatura en x.

El método spline.centroCurv(x)) devuelve las coordenadas {x,y} del centro del círculo de curvatura en x.

Interpolación 2D

Planteamiento

Tenemos una malla en el pano X-Y, de forma que en los nodos tenemos un valor asociado Zi,j = f(xi, yj).
Queremos encontrar una función fi,j(x, y) para ∀ xi ≤ x ≤ xi+1 y ∀ yj ≤ y ≤ yj+1.

Interpolación

Una forma de obtener el valor estimado de f(x,y) es

  1. Calcular los m+1 splines asociados a cada xi, y evaluar f(xi, y)
    2. .
  2. Calcular la spline con los nuevos puntos obtenidos f(xk,y) y calcular sobre ésta el valor estimado de f(x,y).

Pero es procedimiento descrito no nos permite tener una función analítica de interpolación sobre la que podamos, por ejemplo, calcular derivadas parciales o integrar. Ahora el objetivo es buscar una función analítica para cada celda de la malla que se aproxime a la función real. Podemos obtener una expresión de la función de aproximación fi,j(x, y) para ∀ xi ≤ x ≤ xi+1 y ∀ yj ≤ y ≤ yj+1 con el siguiente criterio y procedimiento:

  1. Calculamos todos los splines por filas, para cada fxi(y) aplicado a xi.
    Si sólo nos interesa la celda (i,j) calcularemos los splines para xi y xi+1.
  2. Calculamos los splines por columnas, para cada fyj(x) aplicado a yj. Si sólo nos interesa la celda (i,j) calcularemos los splines para yj e yj+1.
  3. Cada celda tiene definido su contorno por cuatro funciones polinómicas de 3 grado. Nos planteamos buscar una función polinómica de dos variables que contenga en su contorno a estas cuatro funciones, que por simplificar las denotaremos con los subíndices 1, 2, 3 y 4.
    Por tanto, buscamos Z(x,y) tal que cumpla:
    • Z(x=xi, y) = Z1(y) = a1·y3 + b1·y2 + c1·y + d1, donde a1, b1, c1 y d1 los hemos calculado en el paso 1.
    • Z(x=xi+1, y) = Z2(y) = a2·y3 + b2·y2 + c2·y + d2, donde a2, b2, c2 y d2 los hemos calculado en el paso 1.
    • Z(x, y=yj) = Z3(x) = a3·x3 + b3·x2 + c3·x + d3, donde a3, b3, c3 y d3 los hemos calculado en el paso 2.
    • Z(x, y=yj+1) = Z4(x) = a4·x3 + b4·x2 + c4·x + d4, donde a4, b4, c4 y d4 los hemos calculado en el paso 2.

    Esta función es:
    Zi,j(x,y) = s1·x3·y3 + s2·x3·y2 + s3·x2·y3 + s4·x3·y + s5·x2·y2 + s6·x·y3 + s7·x3 + s8·x2·y + s9·x·y5 + s10·y3 + s11·x2 + s12·x·y + s13·y2 + s14·x + s15·y + s16

    Función polinómica con 16 coeficientes que son incógnitas. Este polinomio lo podemos reescribir de dos maneras:
    Zi,j(x,y) = P1(x)·y3 + P2(x)·y2 + P3(x)·y + P4(x) = P5(y)·x3 + P6(y)·x2 + P7(y)·x + P8(x)

    Operando obtememos:
    • P1(x) = s1·x3 + s3·x2 + s6·x + s10
    • P2(x) = s2·x3 + s5·x2 + s9·x + s13
    • P3(x) = s4·x3 + s8·x2 + s12·x + s14
    • P4(x) = s7·x3 + s11·x2 + s14·x + s16
    • P5(y) = s1·y3 + s2·y2 + s4·y + s7
    • P6(y) = s3·y3 + s5·y2 + s8·y + s11
    • P7(y) = s6·y3 + s9·y2 + s12·y + s15
    • P8(y) = s10·y3 + s13·y2 + s15·y + s16
    Condiciones:
    Z1(y) = Zi,j(xi,y) Z2(y) = Zi,j(xi+1,y) Z3(x) = Zi,j(x,yj) Z4(x) = Zi,j(x,yj+1)
    a1 = P1(xi)
    b1 = P2(xi)
    c1 = P3(xi)
    d1 = P4(xi)
    		a2 = P1(xi+1)
    b2 = P2(xi+1)
    c2 = P3(xi+1)
    d2 = P4(xi+1)
    		a3 = P5(yj)
    b3 = P6(yj)
    c3 = P7(yj)
    d3 = P8(yj)
    		a4 = P5(yj+1)
    b4 = P6(yj+1)
    c4 = P7(yj+1)
    d4 = P8(yj+1)
  4. Con todo lo visto, resolvemos el sistema de ecuaciones (los ceros los omitimos):
    xi3 xi2 xi 1
    xi3 xi2 xi 1
    xi3 xi2 xi 1
    xi3 xi2 xi 1
    xi+13 xi+12 xi+1 1
    xi+13 xi+12 xi+1 1
    xi+13 xi+12 xi+1 1
    xi+13 xi+12 xi+1 1
    yj3 yj2 yj 1
    yj3 yj2 yj 1
    yj3 yj2 yj 1
    yj3 yj2 yj 1
    yj+13 yj+12 yj+1 1
    yj+13 yj+12 yj+1 1
    yj+13 yj+12 yj+1 1
    yj+13 yj+12 yj+1 1
    x
    s1
    s2
    s3
    s4
    s5
    s6
    s7
    s8
    s9
    s10
    s11
    s12
    s13
    s14
    s15
    s16
    =
    a1
    b1
    c1
    d1
    a2
    b2
    c2
    d2
    a3
    b3
    c3
    d3
    a4
    b4
    c4
    d4

No ha sido comprobado, puede haber algún error. Es evidente que el conjunto de funciones Zij(x,y) es continua en cualquier dirección. Que es derivable y que la derivada segunda es continua para las líneas xi e yj, pero no está comprobado para el resto de líneas, o para las derivadas parciales 𝜕2f(x,y)/(𝜕x·𝜕y). Pero lo que creo que es claro que la función polinómica con el grado expresado, es la que mejor nos puede encajar para el ajuste de los datos de la malla, siempre que se entienda por "mejor", que la función sea continua, respecto a que lo sean sus derivadas, y 2ª derivadas en cualquier dirección.

Tenemos la forma de evaluar con datos numéricos una función R2🠖R.

Resueltos los coeficientes de la ecuación zi,j(x,y) tendremos una suma de 16 términos de la forma:
Zk=sk·xm·yn.

Podemos calcular:

Derivadas parciales Ver Derivadas parciales

∂Zk/∂x = sk·m·xm-1·yn

∂Zk/∂y = sk·n·xm·yn-1

Gradiente ver documento
rotacional

∇Zk(x,y) = (∂Zk/∂x , ∂Zk/∂y)

Derivadas direccionales rotacional

∇Zk(x,y)·dr = ∂Zk/∂x ·dx + ∂Zk/∂y ·dy

Divergencia y rotacional La divergencia y el rotacional se aplican a campos vectoriales, que se salen del objeto de este documento.
Plano tangente a un punto de una superficie. Vector normal en un punto de una superficie
Integral de volumen de una función 2D en cartesianas Resueltos los coeficientes de la ecuación zi,j(x,y) tendremos una suma de 16 términos Zk=sk·xm·yn. La integral es la suma de las integrales de cada término. La integral de cada término la podemos calcular: \[\int^{y_{j+1}}_{y_j}{\int^{x_{i+1}}_{x_i}{s_k·x^m·y^n·dx}·dy} = \frac{s_k}{(m+1)·(n-1)}·(x_{i+1}^{m+1}-x_i^{m+1})·(y_{j+1}^{n+1}-y_j^{n+1})\]

Integral de superficie de un campo escalar Integrales de superficie
Integrakes de superficie en Wikipedia.
Tratado de integrales de superficie
...
Integral de superficie de un campo vectorial
Integral de línea de un campo escalar Integral de línea en Wikipedia.
Integral de línea de un campo vectorial Integral de línea en Wikipedia.

Resultados:

Tu navegador no admite el elemento <canvas>.