Differenti quadrature per Poisson (1D) (FEM)

Messaggioda feddy » 15/11/2019, 16:37

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.

Immagine

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$.


Immagine


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

Immagine


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?
Avatar utente
feddy
Moderatore
Moderatore
 
Messaggio: 2624 di 5934
Iscritto il: 26/06/2016, 00:25
Località: SISSA

Re: Differenti quadrature per Poisson (1D) (FEM)

Messaggioda feddy » 18/11/2019, 17:04

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
Codice:
Fem1dPoisson


Codice:
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

Immagine


Ultimo bump di feddy effettuato il 18/11/2019, 17:04.
Avatar utente
feddy
Moderatore
Moderatore
 
Messaggio: 2625 di 5934
Iscritto il: 26/06/2016, 00:25
Località: SISSA


Torna a Analisi Numerica e Ricerca Operativa

Chi c’è in linea

Visitano il forum: Nessuno e 1 ospite