19 Approssimazione spettrale e altre applicazioni

Autore/Autrice

Affiliazione

Enrico Bertolazzi

University of Trento, Department of Industrial Engineering

19.1 Approssimazione spettrale

Quanto velocemente decadono i coefficienti

Una funzione sufficientemente regolare può essere approssimata mediante un numero limitato di coefficienti in una base ortogonale. La rapidità del loro decadimento dipende dalla regolarità di f, non dal grado scelto. Se f\in C^m, i coefficienti di Legendre decadono algebricamente come O(k^{-m-1/2}). Se f è analitica e prolungabile all’ellisse di Bernstein E_\rho — l’ellisse di fuochi \pm1 e semiassi di somma \rho>1 — il decadimento è geometrico,

|c_k| \le C\,\rho^{-k}, \tag{19.1}

e l’errore di troncamento al grado N si comporta come \rho^{-N}: è la convergenza spettrale. Una singolarità vicina all’intervallo riduce \rho e quindi rallenta la convergenza senza distruggerla; una discontinuità sull’ intervallo la distrugge, e si ricade nel caso algebrico lento discusso più sotto. Nella pratica i coefficienti si calcolano con quadratura, FFT per Chebyshev oppure ricorrenze e trasformazioni dedicate. ### Una mappa del capitolo

La Figura 19.1 illustra le relazioni tra ortogonalità, approssimazione, quadratura, distribuzione degli zeri e metodi spettrali.

Figura 19.1: Mappa delle relazioni tra famiglie ortogonali, approssimazione, quadratura e metodi spettrali.

Esempio: convergenza dello sviluppo di Legendre di e^x

Calcoliamo i coefficienti della serie di Legendre con una griglia uniforme fine come quadratura composita e misuriamo l’errore massimo su una griglia di validazione. Aumentando il numero di punti di integrazione e di valutazione si può verificare la robustezza della misura. I risultati sono errori osservati, non una maggiorazione analitica.

m<-30000; xq<- -1+(seq_len(m)-0.5)*2/m; wq<-rep(2/m,m)
xx<-seq(-1,1,length.out=2001); truth<-exp(xx); Nmax<-30
Pq<-sapply(0:Nmax,legendre,x=xq)
coef<-colSums(sweep(Pq,1,wq*exp(xq),`*`))*(2*(0:Nmax)+1)/2
errs<-numeric(Nmax+1)
for(n in 0:Nmax){Px<-sapply(0:n,legendre,x=xx)
  vals<-drop(Px%*%coef[1:(n+1)]); errs[n+1]<-max(abs(truth-vals))}
plot(0:Nmax,errs,type="b",pch=16,log="y",lwd=2,col=4,
     xlab="grado massimo N",ylab="max_x |f-p_N|",
     main="Convergenza osservata: espansione di Legendre")

Figura 19.2: Errore massimo misurato per lo sviluppo troncato di Legendre di exp(x).

Il decadimento iniziale rapido è compatibile con l’analiticità di e^x. Quando l’errore si avvicina alla precisione numerica della quadratura e dell’aritmetica, l’andamento può appiattirsi o oscillare.

19.2 I limiti: discontinuità e fenomeno di Gibbs

Un salto produce coefficienti che decadono molto più lentamente e innesca il fenomeno di Gibbs¹. L’aumento del grado restringe la zona oscillante, ma non elimina l’overshoot vicino alla discontinuità. Un filtro sui coefficienti attenua oscillazioni ad alta frequenza, con il costo di allargare la transizione. La Figura 19.3 confronta somme troncate della funzione segno.

m<-40000; xq<- -1+(seq_len(m)-0.5)*2/m; wq<-rep(2/m,m)
f<-sign(xq); xx<-seq(-1,1,length.out=2001); N<-c(12,40)
cols<-c(2,4)
plot(xx,sign(xx),type="s",col="gray50",lty=2,lwd=2,
     xlab="x",ylab="approssimazione",main="Gibbs con polinomi di Legendre",
     ylim=c(-1.35,1.35))
for(j in seq_along(N)){
  n<-N[j]; Pq<-sapply(0:n,legendre,x=xq)
  cf<-colSums(sweep(Pq,1,wq*f,`*`))*(2*(0:n)+1)/2
  Px<-sapply(0:n,legendre,x=xx); lines(xx,drop(Px%*%cf),col=cols[j],lwd=2)
}
legend("bottomright",c("segno(x)",paste0("N=",N)),
       col=c("gray50",cols),lty=c(2,1,1),lwd=2,bty="n")

Figura 19.3: Sviluppo di Legendre della funzione segno: l’aumento del grado restringe ma non annulla le oscillazioni presso il salto.

Un filtro esponenziale smorza i coefficienti con \sigma_k=\exp[-\gamma(k/N)^q] prima di risommare la serie. La scelta dei due parametri non è arbitraria: q è l’ordine del filtro e decide quanto bruscamente si passa dai modi conservati a quelli soppressi, mentre \gamma si prende di norma pari a -\log\varepsilon_{\text{mach}}\approx36 in doppia precisione, così che \sigma_N=\varepsilon_{\text{mach}}: l’ultimo modo viene azzerato esattamente al livello del rumore di arrotondamento, e non oltre. Con q=8 e \gamma=36 si ottiene il filtro usato qui sotto:

m<-40000;xq<- -1+(seq_len(m)-0.5)*2/m;wq<-rep(2/m,m);f<-sign(xq)
xx<-seq(-1,1,length.out=2001);n<-40;Pq<-sapply(0:n,legendre,x=xq)
cf<-colSums(sweep(Pq,1,wq*f,`*`))*(2*(0:n)+1)/2
Px<-sapply(0:n,legendre,x=xx); raw<-drop(Px%*%cf)
sigma<-exp(-36*((0:n)/n)^8); filtered<-drop(Px%*%(cf*sigma))
plot(xx,raw,type="l",col=2,lwd=1.5,xlab="x",ylab="valore",
     main="Gibbs: somma filtrata e non filtrata",ylim=c(-1.35,1.35))
lines(xx,filtered,col=4,lwd=2);lines(xx,sign(xx),col="gray40",lty=2)
legend("bottomright",c("senza filtro","filtro esponenziale","segno(x)"),
       col=c(2,4,"gray40"),lty=c(1,1,2),lwd=c(1.5,2,1),bty="n")

Figura 19.4: Filtro esponenziale: attenua le oscillazioni, rendendo più larga la transizione del salto.

19.3 Metodi spettrali

Collocazione per un problema ai limiti

In un metodo di collocazione spettrale si cerca un polinomio u_N e si impone che il residuo dell’equazione differenziale si annulli in nodi scelti. La matrice di derivazione si costruisce valutando le derivate della base nei nodi; condizioni al bordo e nodi di Chebyshev–Lobatto sono spesso comodi perché includono gli estremi. Per soluzioni lisce, l’errore spettrale può diminuire rapidamente con N, mentre discontinuità e strati limite rallentano la convergenza.

Esempio: per -u''=\pi^2\sin(\pi x), u(-1)=u(1)=0, la soluzione è u(x)=\sin(\pi x). La seguente figura mostra la soluzione calcolata imponendo l’equazione ai nodi interni di Chebyshev–Lobatto. Per ordini alti si usano matrici di derivazione con scalature e formulazioni che controllano il condizionamento.

N<-12; k<-0:N; xn<-cos(pi*k/N)
cvec<-c(2,rep(1,N-1),2)*(-1)^k
DX<-outer(xn,xn,"-"); D<-matrix(0,N+1,N+1)
for(i in 1:(N+1)) for(j in 1:(N+1)) if(i!=j)
  D[i,j]<-(cvec[i]/cvec[j])/DX[i,j]
diag(D)<- -rowSums(D)
D2<-D%*%D; interior<-2:N
rhs<-pi^2*sin(pi*xn[interior])
u<-numeric(N+1)
u[interior]<-solve(-D2[interior,interior],rhs)
xx<-seq(-1,1,length.out=1001); exact<-sin(pi*xx)
lambda<-(-1)^k; lambda[c(1,N+1)]<-lambda[c(1,N+1)]/2
barycentric_eval<-function(z,xn,yn,lambda) {
  vapply(z,function(x) {
    d<-x-xn; hit<-which.min(abs(d))
    if(abs(d[hit])<=8*.Machine$double.eps*max(1,abs(x))) return(yn[hit])
    t<-lambda/d
    sum(t*yn)/sum(t)
  },numeric(1))
}
u_grid<-barycentric_eval(xx,xn,u,lambda)
plot(xx,exact,type="l",lwd=2,col=1,xlab="x",ylab="u(x)",
     main="Collocazione spettrale di -u'' = pi^2 sin(pi x)")
lines(xx,u_grid,lwd=2,lty=2,col=4);points(xn,u,pch=19,col=2,cex=0.6)
legend("topleft",c("soluzione esatta","collocazione","nodi Lobatto"),
       col=c(1,4,2),lty=c(1,2,NA),pch=c(NA,NA,19),bty="n")
cat("errore massimo sulla griglia =",max(abs(u_grid-exact)),"\n")

errore massimo sulla griglia = 2.125982e-07

Figura 19.5: Collocazione Chebyshev–Lobatto della soluzione liscia sin(pi x); indicato l’errore massimo misurato sulla griglia.

19.4 Altre applicazioni

Probabilità, chaos polinomiale e armoniche sferiche

La scelta della misura è centrale anche fuori dall’intervallo finito:

per una variabile uniforme su [-1,1], Legendre è la base naturale;
per una variabile gaussiana standard, i polinomi di Hermite² sono ortogonali rispetto alla densità normale e costituiscono la base del polynomial chaos per propagare incertezza in modelli non lineari;
per coordinate sferiche, i Legendre associati compaiono nelle armoniche sferiche e nella separazione dell’equazione di Laplace³;
i Chebyshev sono usati nell’approssimazione, nell’interpolazione e nei metodi spettrali grazie alla loro relazione con FFT e coseni;
la tridiagonale di Jacobi collega coefficienti, zeri e quadratura, oltre a fornire approssimazioni spettrali di operatori di moltiplicazione.

Gli Hermite e la misura gaussiana

Per la gaussiana standard, la ricorrenza monica di Hermite è He_{n+1}(x)=xHe_n(x)-nHe_{n-1}(x), con \mathbb E[He_n(X)He_m(X)]=n!\,\delta_{nm}. La normalizzazione He_n/\sqrt{n!} rende le basi ortonormali. Ciò consente di espandere una risposta casuale in coefficienti e stimare media e varianza dai coefficienti, senza simulare ogni possibile ingresso.

In sintesi

La famiglia si sceglie in funzione del dominio e della misura: peso uniforme su un intervallo porta a Legendre, peso singolare agli estremi a Jacobi/Chebyshev, densità gaussiana a Hermite. In ogni caso si controllano condizionamento, quadratura usata per i coefficienti, risoluzione delle singolarità e convergenza osservata.

La convergenza spettrale è una proprietà della funzione, non del metodo: 19.1 vale finché f è analitica in un intorno ellittico dell’intervallo.
Una discontinuità produce il fenomeno di Gibbs: l’overshoot non si riduce aumentando N, si restringe soltanto. I filtri lo attenuano al prezzo di allargare la transizione.
Nella collocazione spettrale l’equazione differenziale viene imposta nei nodi; la scelta di nodi di Chebyshev–Lobatto, già incontrati in Capitolo 9, include gli estremi e rende immediate le condizioni al bordo.
Fuori dall’intervallo finito cambia solo la misura: Hermite per la gaussiana, Laguerre per la semiretta. Ricorrenza a tre termini, zeri, quadratura gaussiana valgono identicamente.

19.5 Esercizi

Confrontare gli errori di Legendre e Chebyshev per |x| e per 1/(1+25x^2); per la seconda, stimare \rho in 19.1 a partire dalla posizione dei poli \pm i/5 e confrontare con la pendenza misurata.
Ripetere l’esperimento di Gibbs con una discontinuità spostata in x=0.3 e descrivere l’effetto della base pari.
Per il chaos polinomiale, espandere f(X)=\exp(X) fino al grado 4 nella base di Hermite ortonormale e confrontare media e varianza approssimate con i valori analitici.
Usare nodi di Gauss–Legendre per calcolare \int_{-1}^1 e^{x^2}\,\mathrm dx al crescere di n e separare l’errore di quadratura dall’errore roundoff.