Integración numérica de la función de soporte compacto en un triángulo

10

como sugiere el título, estoy tratando de calcular la integral de una función compacta (polinomio quíntico de Wendland) en un triángulo. Tenga en cuenta que el centro de la función está en algún lugar del espacio tridimensional. Integro esta función en un triángulo arbitrario, pero pequeño ( area<(radius/4)22 ) Actualmente estoy usando la integración descrita por Dunavant, 1985 (p = 19).

Sin embargo, parece que estas reglas de cuadratura no son adecuadas para problemas con soporte compacto. Esto se ve respaldado por el hecho de que cuando integro f(r)=[r1] (una función que es 1 dentro del círculo de radio 1) en un plano que se discretiza usando triángulos, mis resultados (normalizados) están entre 1.001 y 0.897.

Entonces mi pregunta es, ¿existe una regla de cuadratura especializada para este tipo de problema? ¿Funcionaría mejor una regla de integración compuesta de orden inferior?

Desafortunadamente, esta rutina es realmente crítica en mi código, por lo que la precisión es crucial. Por otro lado, necesito hacer esta integración "un par de veces" para un solo paso de tiempo, por lo que el gasto computacional no debería ser demasiado alto. La paralelización no es un problema, ya que ejecutaré la integración en serie.

Gracias de antemano por sus respuestas.

EDITAR: el polinomio quintico de Wendland viene dado por W(q)=[q2]αh3(1q2)4(2q+1)conα=2116π r0R3q=rr0hr0R3

EDIT2: Si es el triángulo bidimensional, entonces quiero calcular con . Entonces en nunca será menor que 0. Tenga en cuenta que la integral es una integral de superficie sobre una superficie 2-D enΔ ω ( r ) d r ω ( r ) = W ( r - r 0ΔΔω(r)drqWR3ω(r)=W(rr0h)qWR3

EDITAR3: Tengo una solución analítica para el problema 1-D (línea). También es posible calcular uno para 2-D (triángulo).

Azrael3000
fuente
¿Podría darnos algunos detalles más de la función que está tratando de integrar? ¿Es solo un polinomio? ¿O un polinomio por partes?
Pedro
Editado según lo solicitado.
Azrael3000

Respuestas:

4

Dado que la función es suave dentro de , pero no de grado fijo (en el plano, es decir), sugeriría usar un esquema adaptativo simple, por ejemplo, la regla trapezoidal con el método de Romberg , en ambas dimensiones.q2

Es decir, si el triángulo es definida por los vértices , y , y usted tiene una rutina de la que se integra a lo largo de la línea de a , se puede hacer lo siguiente (en notación Matlab):y z R 3xyzR3romb(f,a,b)fab

int = romb( @(xi) romb( W , xi , y+(z-y)*(xi-x)./(z-x) ) , x , z );

En romb, no use un número fijo de puntos, pero siga creciendo la tabla hasta que la diferencia entre dos diagonales sucesivas esté por debajo de su tolerancia requerida. Como su función es fluida, esta debería ser una buena estimación de error.

Si partes del triángulo están fuera del dominio de , puede intentar ajustar los límites de integración en el código anterior en consecuencia.W(q)

Puede que esta no sea la forma más eficiente desde el punto de vista computacional de resolver su problema, pero la adaptabilidad le dará mucha más robustez que una regla de grado fijo.

Pedro
fuente
La función es smmoth en todas partes excepto . El vecindario de este punto está causando el problema. q=0
Arnold Neumaier
Ah descomponiéndome en dos problemas 1-D, no es una mala idea en absoluto. Porque hay una cosa que no te he dicho. Tengo una solución analítica en 1-D para poder reemplazar el romb interior por una función analítica. Le daré una oportunidad +1 ya
Azrael3000
@ArnoldNeumaier, lo siento, no veo cómo eso es posible. ¿Podrías explicar?
Pedro
suave como una función de , pero es una función no suave de , y la integración está por encima de , hasta donde yo entiendo la pregunta. La función compuesta es, por lo tanto, una función no suave de . q r r rqqrrr
Arnold Neumaier
1
@Pedro Lo implementé y funciona de maravilla. En realidad, también encontramos una solución analítica hoy. Pero esto es solo para un caso especial que se puede utilizar para reconstruir el general. Eso significa que necesitamos hacer una cierta descomposición del dominio. Dado que el Romberg converge en aproximadamente 4 pasos, creo que debido a esto será más rápido que usar la fórmula analítica. Y según Wikipedia, podemos hacerlo aún mejor que Romberg cuando usamos polinomios racionales. Encontrarás tu nombre en los agradecimientos de mi próximo artículo :) Saludos.
Azrael3000
2

Para una buena visión general de las reglas de cubature, vea "R. Cools, An Encyclopaedia of Cubature Formulas J. Complexity, 19: 445-453, 2003". El uso de una regla fija puede darle la ventaja de que algunas reglas integran polinomios exactamente (como lo hace la cuadratura gaussiana en una dimensión).

Cools es también uno de los principales autores de CUBPACK , un paquete de software para cubicación numérica.

GertVdE
fuente
Creo que el problema aquí es que la función es un polinomio de , pero es una función no lineal en las coordenadas espaciales. La función se suaviza hasta el borde de la función base, pero no es polinómica, excepto a lo largo de los ejes. qqq
Pedro
Esto es correcto Pedro.
Azrael3000
ah ok mi error. lo siento.
GertVdE
2

Las reglas de integración suponen que la función está localmente bien aproximada por un polinomio de bajo grado. Su problema no tiene nada que ver con el soporte compacto. Las funciones de base radial con soporte compacto son suaves en el límite de soporte, y las reglas de cuadratura hasta el orden de suavidad se pueden usar sin problemas. (Las reglas de orden superior no ayudan; por lo tanto, probablemente no debería usar una regla que integre polinomios de grado 5 exactamente).

En su caso, la inexactitud proviene del hecho de que la suposición de buena aproximación polinómica falla en su caso para triángulos cerca de , incluso cuando no contienen .r 0r0r0

q q r r r 0 r rW es suave como una función de , pero es una función no suave de , con un gradiente que se vuelve infinito en el límite . La integración está sobre , y la función compuesta es una función no suave de .qqrrr0rr

Si el triángulo no contiene , la función es pero esto no ayuda ya que la derivada más alta crece muy rápidamente cerca de , y el error de un método de alto orden es proporcional a una derivada de alto orden, por lo tanto muy grande !C i n f r 0r0Cinfr0

El remedio simple es dividir cada triángulo T en un número N_T de sub-triángulos. Puede tomar lejos de y cerca de . Puede averiguar sin conexión qué tan grande debe ser para los triángulos de un diámetro y distancia dados de para alcanzar la precisión deseada. Además, solo debe usar fórmulas de bajo orden cerca de .r 0 N T1 r 0 N T r 0 r 0NT=1r0NT1r0NTr0r0

A medida que se integra sobre un triángulo, pero es tridimensional, el triángulo aparentemente está en .R 3r0R3

Por lo tanto, un remedio más rápido tabularía la integral para como una función de las coordenadas del triángulo (normalizado al rotarlo en un plano bidimensional de modo que un vértice se encuentre en el eje y reflejándolo de tal manera que un segundo vértice se encuentra por encima de él). Esta tabulación debe ser suficientemente detallada para que una interpolación lineal o cuadrática sea lo suficientemente precisa. Pero puede usar el método lento descrito primero para crear esta tabla.x y xr0=0xyx

Otra forma de deshacerse del problema es utilizar una función de base radial con soporte compacto que sea un polinomio en lugar de . Esto es suave en todas partes y fácil de integrar. qq2q

Arnold Neumaier
fuente
Creo que hay un pequeño malentendido. Actualicé la descripción de mi pregunta. De hecho, en la integral nunca puede ser menor que 0. Y no está necesariamente contenida en el triángulo. r 0qr0
Azrael3000
Tu nueva incorporación no tiene sentido para mí. Si entonces también debe ser . ¿O se integra sobre un triángulo 2D en r R 3r0R3rR3r0
R3