Integrando polinômios Lagrange com muitos nós, arredondamento

Dado um conjunto de pontos em , eu gostaria de calcular exatamente. é o polinômio de Lagrange em relação aos pontos com como nó, ou seja, Como esse é um polinômio de grau , eu poderia usar qualquer quadratura gaussiana antiga de grau suficiente. Isso funciona bem se não for muito grande, mas gera resultados com erros de arredondamento para grande . [ - 1 , 1 ] ∫ 1 - 1 L i ( x $\{x_j\}_{j=1}^n$$[-1, 1]$L i x j x i L i ( x ) = ∏ j ≠ i x - x

$$\int_{-1}^{1} L_i(x)\,\text{d} x$$ $L_i$$x_j$$x_i$nn $$L_i(x) = \prod_{j\neq i} \frac{x - x_j}{x_i - x_j}.$$ $n$$n$$n$

Alguma idéia de como evitar isso?

O cálculo de

$$\int_{-1}^{1} L_k(x)\,\text{d} x$$ para os polinômios de Lagrange$L_k$ definidos em uma grade arbitrária$x_k, k=0,\ldots,n$ pode ser realizado pelas duas etapas a seguir:

1. Calcule os pesos de quadratura de Clenshaw-Curtis $w^{\text{cc}}_k$ na grade extrema de Chebyshev $y_k$ para $k=0,\ldots,n$ :

   $$y_k = \cos\left(\frac{k\pi}{n}\right)\\ w^{\text{cc}}_k = \frac{c_k}{n }\Bigg(1-\sum_{j=1}^{\lfloor{n/2}\rfloor} \frac{b_j}{4j^2-1} \cos\left(\frac{2\pi \,j \, k}{n}\right)\Bigg)$$ onde$b_k := 1$para$k=n/2$, caso contrário$b_k:=2$e$c_k:=1$para$k=0$ou$k=n$, caso contrário$c_k:=2$. Veja o artigo de Waldvogel (2006) para mais detalhes.

2. Transformar os pesos $w^{\text{cc}}_k$ à rede arbitrária $x_k, k = 0,\ldots, n$ , através da matriz de transformação $M$ para obter os pesos pretendidos $w_k$ ,

   $$w_k = \sum_{j} M_{kj} w_j^{\text{cc}}$$ onde $$M_{jk} \ = \ L_j(y_k)\,.$$

Em princípio, este é apenas quadratura Clenshaw-Curtis com os valores da função da grelha arbitrária $x_k$ , mas obtidos através de transformação com base (para um refernce geral em Clenshaw-Curtis, ver por exemplo o papel Trefethen).

O algoritmo parece bastante estável, principalmente quando comparado à abordagem de Vandermonde, conforme fornecida na resposta de @Kirill : embora siga as mesmas idéias - gere os pesos de quadratura em uma base conhecida e depois se transforme na nova grade - poderia ter sido esperado, já que a transformação em termos da matriz de Vandermonde é geralmente altamente mal condicionada.

---

Exemplo: Geração de pesos de quadratura de Legendre-Lobatto

Consideramos o exemplo da regra de quadratura de Legendere-Lobatto e comparamos a precisão com a abordagem monomial. Como referência, usamos os pesos em quadratura $w_k^{\text{Leg}}$ obtidos pelo algoritmo Golub-Welsch para diferentes $n$ e calculamos o erro acumulado

$$\epsilon_n = \sum_{k=1}^n \Big(w_k - w_k^{\text{Leg}}\Big)^2$$ Aqui está o resultado: Um observa que os pesos de quadratura de Clenshaw-Curtis são perfeitamente estáveis ​​em toda a faixa considerada de pontos de grade e reproduzem os pesos de Legendre até a precisão da máquina ( $\sqrt{\epsilon} \sim 10^{15}$ ).

---

Exemplo: Geração de fórmulas de quadratura de Newton-Cotes

Consideramos a geração da fórmula de quadratura de Newton-Cotes em grades igualmente espaçadas. Novamente, espera-se um mau condicionamento, pois, em suma, para interpolação polinomial, as grades igualmente espaçadas são baaad.

Na figura a seguir, calculei a soma absoluta dos pesos $\sum_i |w_i|/N$ .

Até, digamos, 50 pontos de grade, o resultado da abordagem monomial e Clenshaw-Curtis concorda. Depois disso, Clenshaw-Curtis se torna melhor - pelo que vale a pena. Uma interpretação direta é que a grade igualmente espaçada arruina tudo para, digamos, $n>10$ . Por volta de $n=50$ , no entanto, a condição da matriz de Vandermonde revida e leva a um resultado ainda pior.

---

Exemplo: quadratura de Guass-Patterson

Este exemplo é devido ao @ NicoSchlömer. Eu não conhecia essas regras até agora, então peguei as abscissas dessa implementação e apliquei a abordagem de Vandermonde e a de Clenshaw-Curtis transformada (onde, como acima, a abordagem de Vandermonde está usando o algoritmo de Björk-Pereyra).

Conforme sugerido no comentário, calculei o erro de integrar uma função constante por

$$\epsilon = \frac{1}{n}\Bigg|2-\sum_{i=1}^n w_i\Bigg|\,,$$ com o seguinte resultado:

A partir dessa imagem, a abordagem transformada de Clenshaw-Curtis parece muito mais eficiente que a abordagem de Vandermonde (pelo menos na aritmética de precisão finita). Ainda assim, Clenshaw-Curtis divide a partir do índice 7, portanto outros métodos devem ser usados.

$b^\top V^{-1}$$b=(2,0,\frac23,0,\frac25,0,\ldots)$

$0\leq x_1<x_2<\cdots <x_n$$V$$b$$0\leq x_1$$O(n^2)$$L_i$$x\mapsto \frac12(x+1)$

$O(\epsilon_{\mathrm{mach}})$

module VandermondeInverse

using SpecialMatrices

function main(n=8)
  X = Rational{BigInt}[k//(n-1) for k=0:n-1]
  # X = convert(Vector{Rational{BigInt}}, linspace(-1, 1, n))
  x = convert(Vector{Float64}, X)

  A = convert(Matrix{Rational{BigInt}}, Vandermonde(X))
  b = [i%2==0 ? 2//(i+1) : 0 for i=0:n-1]
  println("Norm: ", norm(A, Inf))
  println("Norm of inverse: ", norm(inv(A), Inf))
  println("Condition number: ", cond(convert(Matrix{Float64}, A)))
  ans = A'\b
  println("True answer: ", ans)

  B = convert(Matrix{Float64}, A)
  c = convert(Vector{Float64}, b)

  println("Linear solve: ", norm((B'\c - ans)./ans, Inf))

  d = vec(c')
  for k=1:n, l=n:-1:k+1
    d[l] -= x[k]*d[l-1]
  end

  for k=n-1:-1:1, l=k:n
    if l > k
      d[l] /= x[l]-x[l-k]
    end
    if l < n
      d[l] -= d[l+1]/(x[l+1] - x[l-k+1])
    end
  end
  println("Neville elimination: ", norm((d-ans)./ans, Inf))

  nothing
end

end

V = VandermondeInverse

Resultado:

julia> V.main(14)
Norm: 14.0
Norm of inverse: 1.4285962612120493e10
Condition number: 5.2214922998851654e10
True answer: Rational{Int64}[3202439130233//2916000,-688553801328731//52390800,19139253128382829//261954000,-196146528919726853//785862000,6800579086408939//11642400,-43149880138884259//43659000,32567483200938127//26195400,-7339312362348889//6237000,48767438804485271//58212000,-69618881108680969//157172400,44275410625421677//261954000,-2308743351566483//52390800,11057243346333379//1571724000,-209920276397//404250]
Linear solve: 1.5714609387747318e-8
Neville elimination: 1.3238218572356314e-15

Se Xnão for positivo como neste teste, parece que os erros relativos são da mesma ordem que com uma resolução linear regular.

$b^\top V^{-1}$$L_i$$L_i(x_j)=\delta_{ij}$$\alpha_{jk}$$L_k$

$$L_k(x) = \sum_{j,k} \alpha_{j,k}x^j = (1,x,x^2,\ldots,x^n)^\top (\alpha_{0k},\ldots,\alpha_{nk}),$$ $L$ $$L = \begin{pmatrix} \alpha_{00}& \cdots & \alpha_{0n}\\\vdots &&\vdots\\ \alpha_{n0}&\cdots&\alpha_{nn} \end{pmatrix}.$$ $L_k$$L$$(1,x,\ldots,x^n)$ $$(1,x,x^2,\ldots,x^n)L = (L_0(x),L_1(x),\ldots,L_n(x)).$$ $L_k(x_j)=\delta_{jk}$ $$\begin{pmatrix} 1&x_0&x_0^2&\cdots&x_0^n\\ \vdots\\ 1&x_n&x_n^2&\cdots&x_n^n \end{pmatrix} L = I,$$ $L = V^{-1}$$V$$x_j$

$\int_{-1}^{1}x^k\,\mathrm{d}x = \frac{1+(-1)^k}{k+1}$

$$\int_{-1}^{1}L_k(x)\,\mathrm{d}x = \sum_j \alpha_{jk}\frac{1+(-1)^k}{k+1} = (2,0,\tfrac23,0,\tfrac25,0,\ldots)^\top (\alpha_{0k},\ldots,\alpha_{nk}).$$ $n+1$$k=0\ldots n$$(2,0,\tfrac23,0,\ldots)^\top L$$L=V^{-1}$

Calcule os produtos dos nominadores e denominadores primeiro e depois divida uma vez. Os dois produtos devem ser da mesma ordem de grandeza, para que não haja erros de arredondamento significativos. Além disso, você obtém o benefício adicional de maior velocidade, devido ao número reduzido de cálculos de ponto flutuante.