Integração numérica da função compactamente suportada em um triângulo

10

como o título sugere, estou tentando calcular a integral de uma função compactamente suportada (polinômio quintic de Wendland) em um triângulo. Observe que o centro da função está em algum lugar no espaço 3D. Integro essa função em um triângulo arbitrário, mas pequeno ( ) Atualmente, estou usando a integração descrita por Dunavant, 1985 (p = 19).area<(radius/4)22

Parece, no entanto, que essas regras de quadratura não são adequadas para problemas com suporte compacto. Isso é suportado pelo fato de que, quando integro (uma função que é 1 dentro do círculo do raio 1) em um plano discretizado usando triângulos, meus resultados (normalizados) ficam entre 1.001 e 0.897.f(r)=[r1]

Então, minha pergunta é: existe uma regra de quadratura especializada para esse tipo de problema? Uma regra de integração composta de ordem inferior funcionaria melhor?

Infelizmente, essa rotina é realmente crítica no meu código, portanto a precisão é crucial. Por outro lado, preciso fazer essa integração "algumas vezes" por uma única etapa, para que as despesas computacionais não sejam muito altas. Paralelização não é um problema, pois executarei a própria integração em série.

Agradecemos antecipadamente por suas respostas.

EDIT: O polinômio quintic de Wendland é dado por comα=21W(q)=[q2]αh3(1q2)4(2q+1)α=2116π r0R3q=rr0hr0R3

EDIT2: Se é o triângulo bidimensional, desejo calcular com . Portanto, em nunca será menor que 0. Observe que a integral é uma integral de superfície sobre uma superfície 2-D emô ω ( r ) d r ω ( r ) = W ( r - r 0ΔΔω(r)drqWR3ω(r)=W(rr0h)qWR3

EDIT3: Eu tenho uma solução analítica para o problema 1-D (linha). Também é possível calcular um para 2-D (triângulo).

Azrael3000
fonte
Você poderia nos dar mais alguns detalhes da função que você está tentando integrar? É apenas um polinômio? Ou um polinômio por partes?
Pedro
Editado conforme solicitado.
Azrael3000

Respostas:

4

Como a função é suave dentro de , mas não de grau fixo (no plano, isto é), sugiro usar um esquema adaptativo simples, por exemplo, a Regra Trapezoidal com o método de Romberg , em ambas as dimensões.q2

Ou seja, se o triângulo é definido pelo vértices , e , e você tem uma rotina que integra ao longo da linha a partir de , você pode fazer o seguinte (em notação Matlab):y z R 3xyzR3romb(f,a,b)fab

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

Em romb, não use um número fixo de pontos, mas continue aumentando a tabela até que a diferença entre duas diagonais sucessivas esteja abaixo da tolerância necessária. Como sua função é suave, essa deve ser uma boa estimativa de erro.

Se partes do triângulo estão fora do domínio de , você pode tentar ajustar os limites de integração no código acima de acordo.W(q)

Essa pode não ser a maneira mais eficiente em termos computacionais de resolver seu problema, mas a adaptabilidade fornecerá muito mais robustez do que uma regra de grau fixo.

Pedro
fonte
A função é smmoth em todos os lugares, exceto . A vizinhança deste ponto está causando o problema. q=0
Arnold Neumaier 22/03
Ah, decompondo-me em dois problemas 1-D, não é uma má idéia. Porque há uma coisa que eu não te disse. Eu tenho uma solução analítica em 1-D para poder substituir o romb interior por uma função analítica. Já vou dar uma chance a isso com +1
Azrael3000
@ArnoldNeumaier, desculpe, não vejo como isso é possível. Você poderia explicar?
Pedro
suave como uma função de , mas é uma função não suave de , e a integração é superior a , tanto quanto eu entendi a pergunta. A função composta é, portanto, uma função não suave de . q r r rqqrrr
Arnold Neumaier 22/03/12
11
@Pedro eu o implementei e funciona como um encanto. Na verdade, também encontramos uma solução analítica hoje. Mas isso é apenas para um caso especial que pode ser usado para reconstruir o geral. Isso significa que precisamos fazer alguma decomposição de domínio. Como o Romberg converge em cerca de 4 etapas, acho que, por isso, será mais rápido do que usar a fórmula analítica. E, de acordo com a Wikipedia, ainda podemos fazer melhor do que Romberg ao usar polinômios racionais. Você encontrará seu nome nos agradecimentos do meu próximo artigo :) Felicidades.
Azrael3000
2

Para uma boa visão geral das regras de cubatura, consulte "R. Cools, Uma Enciclopédia de Fórmulas de Cubatura J. Complexity, 19: 445-453, 2003". Usando uma regra fixa, você pode obter a vantagem de que algumas regras integram polinômios exatamente (como a quadratura gaussiana faz em uma dimensão).

Cools também é um dos principais autores do CUBPACK , um pacote de software para cubatura numérica.

GertVdE
fonte
Eu acho que o problema aqui é que a função é um polinômio de , mas é uma função não linear nas coordenadas espaciais. A função é suavizada até a borda da função base, mas não polinomial, exceto ao longo dos eixos. qqq
Pedro
Está correto, Pedro.
Azrael3000
Ah ok. meu erro. desculpa.
GertVdE
2

As regras de integração assumem que a função é localmente bem aproximada por um polinômio de baixo grau. Seu problema não tem nada a ver com suporte compacto. As funções de base radial compactamente suportadas são suaves no limite do suporte, e regras de quadratura até a ordem da suavidade podem ser usadas sem problemas. (Regras de ordem superior não ajudam; portanto, você provavelmente não deve usar uma regra que integre exatamente os polinômios do grau 5).

No seu caso, a imprecisão deriva do fato de que a suposição de boa aproximação polinomial falha no seu caso para triângulos próximos a , mesmo quando eles não contêm .r 0r0r0

q q r r r 0 r rW é suave como uma função de , mas é uma função não suave de , com um gradiente que se torna infinito no limite . A integração ultrapassa , e a função composta é uma função não suave de .qqrrr0rr

Se o triângulo não contém , a função é mas isso não ajuda, pois a derivada mais alta cresce muito rapidamente perto de , e o erro de um método de alta ordem é proporcional a uma derivada de alta ordem, portanto, muito grande !C i n f r 0r0Cinfr0

O remédio simples é dividir cada triângulo T em um número N_T de sub-triângulos. Você pode levar para longe de e para perto de . Você pode descobrir off-line qual o tamanho de para triângulos de um determinado diâmetro e distância de para alcançar a precisão desejada. Além disso, você deve usar apenas fórmulas de ordem baixa perto de .r 0 N T » 1 r 0 N T R 0 R 0NT=1r0NT1r0NTr0r0

À medida que você integra sobre um triângulo, mas é tridimensional, aparentemente o triângulo está em .R 3r0R3

Um remédio mais rápido, portanto, tabulará a integral para em função das coordenadas do triângulo (normalizada girando-a em um plano bidimensional, de modo que um vértice esteja no eixo e refletindo-a de modo que um segundo está acima dele). Essa tabulação deve ser suficientemente detalhada para tornar uma interpolação linear ou quadrática precisa o suficiente. Mas você pode usar o método lento descrito primeiro para criar esta tabela.x y xr0=0xyx

Outra maneira de se livrar do problema é usar uma função de base radial compactamente suportada que é um polinômio em vez de . Isso é suave em qualquer lugar e fácil de integrar. qq2q

Arnold Neumaier
fonte
Eu acho que há um pequeno mal-entendido. Eu atualizei a descrição da minha pergunta. De fato, na integral nunca pode ser menor que 0. E não está necessariamente contido no triângulo. r 0qr0
Azrael3000
Sua nova adição não faz sentido para mim. Se , então deve ser . Ou você integra sobre um triângulo 2D em ? - Não assumi que esteja no triângulo. Acabei de adicionar mais um momento à minha resposta. r R 3 r 0r0R3rR3r0
Arnold Neumaier 22/03/12
Sim, é correto que eu integre um triângulo 2D em . R3
Azrael3000