Integrando polinomios de Lagrange con muchos nodos, redondeo

10

Dado un conjunto de puntos en [ - 1 , 1 ] , me gustaría calcular 1 - 1 L i ( x ){xj}j=1n[1,1] exactamente. L i es el polinomio de Lagrange con respecto a los puntos x j con x i como nodo, es decir, L i ( x ) = j i x - x j

11Li(x)dx
Lixjxi Como este es un polinomio de gradon, podría usar cualquier cuadratura gaussiana de grado suficiente. Esto funciona bien sinno es demasiado grande, pero conduce a resultados defectuosos por errores de redondeo parangrande.
Li(x)=jixxjxixj.
nnn

¿Alguna idea de cómo evitarlos?

Nico Schlömer
fuente
3
Esto depende de dónde estén las 's, pero ¿ha verificado que sus L i se comportan bien? En el peor de los casos, con x j distribuido uniformemente, obtienes el fenómeno Runge ( L i es oscilatorio y grande), en cuyo caso no son realmente errores de redondeo los que causan problemas. xjLixjLi
Kirill
2
Además, nitpick: dividir entre números pequeños es una operación bien condicionada, es más bien la substracción subsiguiente de números grandes casi iguales que está mal condicionada y conduce a la inestabilidad numérica.
Kirill
Parece que estás tratando de calcular dondeVes la matriz de Vandermonde dexj's. ¿Puedes decir cuál es el número de condición deV? (2,0,23,0,25,0,)V1VxjV
Kirill

Respuestas:

5

El cálculo de

11Lk(x)dx
para los polinomios de LagrangeLk definidos en una cuadrícula arbitrariaxk,k=0,,n puede realizarse mediante los dos pasos siguientes:

  1. Calcule los pesos de la cuadratura de Clenshaw-Curtis wkcc en la cuadrícula extrema de Chebyshev yk para k=0,,n :

    yk=cos(kπn)wkcc=ckn(1j=1n/2bj4j21cos(2πjkn))
    dondebk:=1parak=n/2, de lo contrariobk:=2, yck:=1parak=0ok=n, de lo contrariock:=2. Ver el artículo de Waldvogel (2006) para más detalles.

  2. Transforme los pesos wkcc en la cuadrícula arbitraria xk,k=0,,n , a través de la matriz de transformación M para obtener los pesos buscados wk ,

    wk=jMkjwjcc
    donde
    Mjk = Lj(yk).

xk

El algoritmo parece ser bastante estable, particularmente cuando se lo compara con el enfoque de Vandermonde según lo provisto en la respuesta de @Kirill : aunque sigue las mismas ideas: generar los pesos en cuadratura de manera conocida y luego transformarlos en la nueva cuadrícula. podría esperarse ya que la transformación en términos de la matriz de Vandermonde suele estar muy mal condicionada.


Ejemplo: generación de pesos en cuadratura Legendre-Lobatto

wkLegn

ϵn=k=1n(wkwkLeg)2
ingrese la descripción de la imagen aquíϵ1015


Ejemplo: generación de fórmulas de cuadratura de Newton-Cotes

Consideramos la generación de la fórmula de la cuadratura de Newton-Cotes en cuadrículas equidistantes. Una vez más, uno espera un mal acondicionamiento, ya que, en resumen, para la interpolación polinómica, las redes espaciadas equitativamente son malas.

i|wi|/Ningrese la descripción de la imagen aquí

n>10n=50


Ejemplo: cuadratura de Guass-Patterson

Este ejemplo se debe a @ NicoSchlömer. No conocía estas reglas hasta ahora, así que tomé las abscisas de esta implementación y apliqué tanto el enfoque de Vandermonde como el enfoque transformado de Clenshaw-Curtis (donde, como anteriormente, el enfoque de Vandermonde está utilizando el algoritmo Björk-Pereyra).

ϵ=1n|2i=1nwi|,

ingrese la descripción de la imagen aquí A partir de esta imagen, el enfoque transformado de Clenshaw-Curtis parece mucho más eficiente que el enfoque de Vandermonde (al menos en aritmética de precisión finita). Aún así, Clenshaw-Curtis se desglosa a partir del índice 7, por lo que se deben utilizar otros métodos.

davidhigh
fuente
wk
Olvidé mencionar que esto es con los puntos de Gauss-Patterson, github.com/nschloe/quadpy/blob/master/quadpy/line_segment/… .
Nico Schlömer
@ NicoSchlömer: ok, eso es interesante, porque a excepción de la aritmética de precisión arbitraria (que la implementación de julia parece ofrecer), es difícil razonar cómo el enfoque de Vandermonde debería ser más preciso (ya que Vandermonde es uno de los mejores ejemplos para un matriz mal acondicionada). Haré una edición de la pregunta y consideraré la cuadratura de Gauss-Patterson.
davidhigh
Gracias David por investigar todo esto. ¡Es realmente interesante! Creo que valdría la pena agregar esta publicación para incluir una comparación con la cuadratura de Gauss-Legendre ordinaria del grado apropiado. Supongo que funciona como su enfoque Clenshaw-Curtis. Además, poner el código en github más o menos y vincularlo desde aquí será útil para cualquiera que esté investigando esto en el futuro. Si vincula la respuesta mejor votada, haré que esta sea la "correcta" debido a las interesantes ideas.
Nico Schlömer
@ NicoSchlömer: gracias. ¿Qué quieres decir en comparación con Gauss-Legendre? (porque los pesos Gauss-Legendre se reproducen con precisión de máquina). La comparación entre Leg. y Trefethen realizó CC, con el resultado de que la precisión de CC es a menudo comparable. Lo interesante sería estudiar el rendimiento de diferentes cuadrículas personalizadas y compararlo con Legendre o Clenshaw-Curtis. Además, vinculé la respuesta más votada.
davidhigh
10

bV1b=(2,0,23,0,25,0,)

0x1<x2<<xnVb0x1O(n2)Lix12(x+1)

O(ϵ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

Salida:

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

Si Xno es positivo como en esta prueba, parece que los errores relativos son del mismo orden que con una resolución lineal regular.

bV1LiLi(xj)=δijαjkLk

Lk(x)=j,kαj,kxj=(1,x,x2,,xn)(α0k,,αnk),
L
L=(α00α0nαn0αnn).
LkL(1,x,,xn)
(1,x,x2,,xn)L=(L0(x),L1(x),,Ln(x)).
Lk(xj)=δjk
(1x0x02x0n1xnxn2xnn)L=I,
L=V1Vxj

11xkdx=1+(1)kk+1

11Lk(x)dx=jαjk1+(1)kk+1=(2,0,23,0,25,0,)(α0k,,αnk).
n+1k=0n(2,0,23,0,)LL=V1
Kirill
fuente
Gracias por la elaborada respuesta! Para agregar aún más claridad, ¿podría dar más detalles sobre cómo llevar el problema a la forma de Vandermonde?
Nico Schlömer
@ NicoSchlömer Claro, ver edición. Gracias por la pregunta, no tenía idea de que este algoritmo existiera.
Kirill
n=16n=31
@ NicoSchlömer Hm, no puedo reproducir eso: con el código anterior, V.main(32)produce una respuesta sensata en aproximadamente un segundo en mi computadora portátil (mientras uso solo un poco de memoria). Los números ni siquiera son tan grandes, el numerador más grande tiene 54 dígitos, por lo que sospecho que algo más te va mal. ¿Puedes publicar una esencia, porque tengo curiosidad por ver cómo falla?
Kirill
1
@ NicoSchlömer Si te refieres a la forma en que comienza a usar BigFloat al imprimir los errores relativos para la salida, no estoy muy seguro de cuál es la regla. Pero me aseguré de que use Float64para d: consultar con @show typeof(d). Avísame si encuentras más problemas con él.
Kirill
2

Calcule primero los productos de los nominadores y los denominadores y luego divídalos una vez. Los dos productos deben ser del mismo orden de magnitud, por lo que no debe haber errores significativos de redondeo. También obtiene el beneficio adicional de una mayor velocidad, debido a la reducción del número de cálculos de coma flotante.

borrar
fuente