Differenti quadrature per Poisson (1D) (FEM)

[feddy]

Ciao a tutti,

questo post era partito come una domanda, poi però scrivendola credo di essermi risposto da solo. Mi farebbe piacere sapere se questa mia "riflessione" sia corretta.


Considerando l'equazione di Poisson $-u''=f$ con condizioni al bordo essenziali, sappiamo che la formulazione debole,una volta ristretti ad un sottospazio $V_h$ finito dimensionale di $H_0^1$ con base data dalle hat-functions $\phi$

trovare $u \in V_h$ t.c. $- \sum_{i=1}^{M-1} u_i \langle \partial_x \phi_i, \partial_x \phi_j \rangle =\langle f, \phi_j \rangle$ valga per ogni $j$

Venendo al mio problema, ossia la quadratura del termine di destra. Prendendo una griglia non equispaziata e come termine di destra $f(x)=cos(2*pi*x)$, ritrovo esattamente la soluzione analitica usando una formula di quadratura gaussiana e praticamente l'errore puntuale rispetto alla soluzione analitica è nullo in aritmetica di macchina, come è lecito aspettarsi, visto che questa è noto per l'equazione di Poisson. Ovviamente, nel caso equispaziato il risultato in termini di accuratezza è identico.

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

Considerando ora $f(x)=e^{-10^4 (x-0.5)^2}$, come si può vedere questa ha un comportamente molto simile ad una delta di Dirac centrata a $x=\frac{1}{2}$.

Utilizzando una griglia equispaziata e formula di quadratura gaussiana, ottengo ancora la soluzione analitica. Lo stesso accade per una griglia non equispaziata. Tutto questo senza prendere "tanti" punti.

- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

Supponiamo ora però di utilizzare la formula del trapezio per calcolare ogni $\langle f, \phi_i \rangle$.In questo caso da quello che ho intuito a lezione, una griglia equispaziata o meno fa la differenza

Nel nodo i-esimo si ha che $f_i= \langle f, \phi_i \rangle = \int_{x_{i-1}}^{x_{i+1}} f(x) \phi \approx f(x_i) \frac{h_{i-1}+h_i}{2}$.

Ora, se la griglia fosse equispaziata, e diciamo non troppo fine, allora chiaramente non mi consentirà di catturare il salto che fa la funzione e perciò la valutazione di $f$ nei nodi sarà praticamente nulla.

Con 41 punti, questo è il risultato.



Incrementando il numero di punti riesco ad avvicinarmi sempre di più alla soluzione analitica. Il seguente grafico di convergenza (loglog), (in $x$ il numero di punti, in $y$ la norma infinito dell'errore), mostra che prendendo punti a sufficienza l'ordine è $2$.





Se invece infittisco la griglia intorno a dove ho il picco, ossia a $x= \frac{1}{2}$, ottengo già un risultato apprezzabile. Infatti, sempre prendendo 41 punti come nella prima delle due figure qui sopra, si ottiene




Quindi, in sostanza, quando uso la routine quad dalla libreria scipy questa fa già il lavoro di suddivisione e pertanto nonostante la griglia non sia troppo fine, calcola comunque molto bene l'integrale.


C'è altro che mi sfugge?

[feddy]

Nel frattempo, ho trovato in rete un modo per adattare la griglia grazie ad una stima a posteriori (un po' come si adatta il passo nelle ODE). Sicuramente sarà la versione più naive per farlo, ma per ora mi accontento di questo.

Riporto per completezza solo il loop (qui Fem1dPoisson è la routine che calcola l'approssimazione di $u$ risolvendo il sistema lineare nella discussione sopra su una generica griglia nonequispaziata).

Uscito dal ciclo calcolo quindi la soluzione tramite la routine

Fem1dPoisson

M = 10
x = np.linspace(0,1,M+1)
while (M<25):
    
    rho = np.zeros(M)
    for k in range(0,M):
        h = x[k+1] - x[k]
        a = rhs(x[k])
        b = rhs(x[k+1])
        t = 0.5*h*(a**2 + b**2)
        rho[k] = t*h**2
        
    alfa = 0.9
    for i in range(len(rho)):
        if rho[i]>alfa*np.max(rho):
            x = np.append(x,(x[i+1]+x[i])/2)
    
    x = np.sort(x)
    M = len(x) - 1
    

u_int = Fem1dPoisson(M,x,rhs,quadr)

E questo è il risultato