Risoluzione Numerica di PDE

Sia \Omega \subseteq \mathbb{R}^d un aperto connesso, u: \overline{\Omega} \to \mathbb{R} e L[u] un operatore differenziale. Consideriamo il problema:

L[u] = f

con f: \overline{\Omega} \to \mathbb{R} nota e sufficientemente regolare.

In questo capitolo ci concentreremo su PDE del secondo ordine in due variabili spaziali (x,y) \in \Omega \subseteq \mathbb{R}^2:

L[u] = A \frac{\partial^2 u}{\partial x^2} + B \frac{\partial^2 u}{\partial x \partial y} + C \frac{\partial^2 u}{\partial y^2} + D \frac{\partial u}{\partial x} + E \frac{\partial u}{\partial y} + Fu

con A, B, C, D, E, F: \Omega \to \mathbb{R}.

Le condizioni al bordo fissano in qualche modo il valore di u sulla frontiera di \Omega.

Possiamo considerare una conica associata al problema:

c(x,y) \coloneqq Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0

\mathcal{C} := \{ (x,y) \in \mathbb{R}^2 \mid c(x,y) = 0 \}

\begin{vmatrix} a & b/2 & c/2 \\ b/2 & d & e/2 \\ c/2 & e/2 & f \end{vmatrix} \neq 0

Tipi di problemi

Problemi ellittici

Consideriamo l’operatore:

\Delta u := \operatorname{tr}(\nabla^2 u) = \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}

detto Laplaciano. In particolare l’equazione differenziale \Delta u = 0 è detta equazione di Laplace.

\begin{cases}
    \dfrac{\partial u}{\partial x} = \dfrac{\partial v}{\partial y} \\[1em]
    \dfrac{\partial u}{\partial y} = -\dfrac{\partial v}{\partial x}
\end{cases}

\begin{cases}
    \dfrac{\partial^2 u}{\partial x^2} = \dfrac{\partial^2 v}{\partial x \partial y} \\[1em]
    \dfrac{\partial^2 u}{\partial y^2} = -\dfrac{\partial^2 v}{\partial y \partial x}
\end{cases} \implies \Delta u = 0 \text{ e } \Delta v = 0

Equazione di Poisson

\Delta u = f(x,y)

A = \iint_{\Omega} \sqrt{1 + \left( \frac{\partial u}{\partial x} \right)^2 + \left( \frac{\partial u}{\partial y} \right)^2} \, \mathrm dx \mathrm dy

Problemi parabolici

Un esempio classico è dato dai problemi dipendenti dal tempo, come l’equazione del calore:

\gamma \frac{\partial u(x,t)}{\partial t} - \frac{\partial^2 u(x,t)}{\partial x^2} = 0, \quad 0 < x < \ell, t \ge 0

con condizioni iniziali e al bordo:

\begin{cases} u(x,0) = g(x) & x \in [0,\ell] \\ u(0,t) = L \\ u(\ell,t) = R \end{cases}

In più dimensioni, l’equazione assume la forma:

\gamma \frac{\partial u(\mathbf{x},t)}{\partial t} - \Delta_{\mathbf{x}} u = 0

Problemi iperbolici

Questi problemi sono spesso legati alle “leggi di conservazione”.

\begin{cases} \dfrac{\partial^2 u}{\partial t^2} - \lambda \dfrac{\partial^2 u}{\partial x^2} = 0 & \lambda > 0 \\[1em] u(\cdot, 0) = u_0, \quad \dfrac{\partial u}{\partial t} (\cdot, 0) = v_0 \\[1em] u(0,t) = u(1,t) = 0 \end{cases}

\nabla^2 v = -\frac{\omega^2}{\lambda} v

\begin{cases} \dfrac{\partial u}{\partial t} - c \dfrac{\partial u}{\partial x} = f(x,t) \\[1em] u(\cdot, 0) = g \end{cases}

Discretizzazioni

L’obiettivo è approssimare la soluzione u di una PDE tramite un sistema algebrico.

Introduzione

Sia n > 0. Approssimiamo u(x) su una griglia con n punti, vorremmo che u_i \approx u(x_i). Questo trasforma l’operatore differenziale L in un operatore lineare L_n di dimensione n:

L u = f \implies L_n \mathbf{u} = \mathbf{f}

In questo ambito ci interessa studiare:

• Come costruire L_n partendo da L.

• L’ordine di convergenza.

• Quali proprietà di L vengono ereditate da L_n.

Metodo di Galerkin

Supponiamo che l’operatore L sia definito su uno spazio di Hilbert H. Il problema L[u] = f può essere formulato in forma debole: cerchiamo u \in H tale che

\langle L[u] - f, v \rangle = 0 \quad \forall v \in H

Sia \{ \phi_i \}_{i \ge 0} una base di H. Esprimiamo u come combinazione lineare degli elementi della base: u = \sum_{i \ge 0} u_i \phi_i. Sostituendo nella formulazione debole e scegliendo v = \phi_j si ottiene:

\langle L \left[ \sum_{i \ge 0} u_i \phi_i \right] - f, \phi_j \rangle = 0 \quad \forall j

che si semplifica nel sistema lineare:

\sum_{i \ge 0} u_i \langle L[\phi_i], \phi_j \rangle = \langle f, \phi_j \rangle \quad \forall j

Problema modello: Laplaciano 1D

Consideriamo il seguente problema modello:

\begin{cases} u''(x) = f(x) & x \in (0,1) \\ u(0) = \alpha \\ u(1) = \beta \end{cases}

Per il teorema fondamentale del calcolo integrale, possiamo scrivere:

u(t) = u(0) + u'(0) t + \int_0^t \int_0^s f(x) \, \mathrm dx \mathrm ds

Questo approccio richiede un’elevata regolarità per le funzioni coinvolte, ad esempio se assumiamo che f \in C^2 allora u \in C^4.

Discretizzazione

Consideriamo un problema della forma L[u] = f. Sia u(x_i) \approx u_i. Se l’operatore differenziale è L[u] = u'', per Taylor abbiamo:

\begin{align}
u''(x) &= \frac{u(x-h) - 2u(x) + u(x+h)}{h^2} + \frac{h^2}{24}(u^{(4)}(\xi) + u^{(4)}(\eta)) \\
&= \frac{u(x-h) - 2u(x) + u(x+h)}{h^2} + h^2 \tau
\end{align}

con |\tau| \le \frac{1}{12} M_4, dove M_4 = \sup_x |u^{(4)}(x)|.

Dunque, se consideriamo la discretizzazione L[u]=u'' \rightsquigarrow L_n u = \frac{1}{h^2}(u_{i-1} - 2u_i + u_{i+1})_i, da cui otteniamo:

L_n u = f - h^2 \tau

Definito in questo modo, L_n: \mathbb{R}^{n+2} \to \mathbb{R}^n è una matrice n \times (n+2):

L_n = \frac{1}{h^2} \begin{bmatrix}
  1 & -2 & 1 & & & \\
  & 1 & -2 & 1 & & \\
  & & \ddots & \ddots & \ddots & \\
  & & & 1 & -2 & 1
\end{bmatrix}

Il vettore di componenti h^2 \tau = (L[u]|_{x=x_i} - L_n u)_i è detto errore locale di discretizzazione. In generale, l’approssimazione si dice consistente di ordine p se \exists \gamma > 0 tale che:

|L[u]|_{x_i} - L_n u| \le \gamma h^p

Possiamo togliere la prima e l’ultima colonna e riscrivere il sistema come segue con una matrice quadrata A_n:

-\frac{1}{h^2} \begin{bmatrix}
  -2 & 1 & & \\
  1 & -2 & 1 & \\
  & \ddots & \ddots & \ddots \\
  & & 1 & -2 & 1
\end{bmatrix} \mathbf{u} = \mathbf{f} - h^2 \mathbf{\tau} - \frac{1}{h^2} \begin{bmatrix}
  a \\
  0 \\
  \vdots \\
  0 \\
  b
\end{bmatrix}

e chiamiamo A_n := -\frac{1}{h^2} \begin{bmatrix} \dots \end{bmatrix} la nuova matrice discretizzata. A questo punto possiamo considerare un nuovo problema associato al precedente senza l’errore locale:

A_n \mathbf{v} = \mathbf{f} - \frac{1}{h^2}(a\mathbf{e}_1 + b\mathbf{e}_n) \quad (1)

Nell’incognita \mathbf{v}, a questo punto introduciamo \mathbf{\epsilon}^{(n)} = \mathbf{v} - \mathbf{u} come l’errore globale, che ci permette di quantificare l’accuratezza del metodo. Vorremmo che \|\mathbf{\epsilon}^{(n)}\|_\infty \to 0 in modo che infittendo la discretizzazione possiamo sperare di convergere alla soluzione esatta.

Per stimare \|\mathbf{\epsilon}^{(n)}\|_\infty sottraiamo (1) e (2), dove (2) è:

A_n \mathbf{u} = \mathbf{f} - h^2 \mathbf{\tau}^{(n)} - \frac{1}{h^2}(a\mathbf{e}_1 + b\mathbf{e}_n) \quad (2)

Otteniamo:

A_n \mathbf{\epsilon}^{(n)} = h^2 \mathbf{\tau}^{(n)} \implies \mathbf{\epsilon}^{(n)} = h^2 A_n^{-1} \mathbf{\tau}^{(n)}

\implies \|\mathbf{\epsilon}^{(n)}\| \le h^2 \|A_n^{-1}\| \|\mathbf{\tau}^{(n)}\|

Studiamo ora la convergenza in norma 2 e in norma \infty.

Analisi errore in norma 2

Vorremmo trovare gli autovalori della matrice A_n. Osserviamo che A_n = 2I - H_n con H_n = \text{trid}(1,0,1). Gli autovalori di A_n sono:

\lambda_j(A_n) = 2 - 2\cos\left(\frac{j\pi}{n+1}\right) = 4 \sin^2\left(\frac{j\pi}{2(n+1)}\right)

per j=1, \dots, n. L’inversa ha autovalori \lambda_j(A_n^{-1}) = 1/\lambda_j(A_n). La norma 2 di una matrice simmetrica è il suo raggio spettrale.

\|A_n^{-1}\|_2 = \max_j \frac{1}{\lambda_j(A_n)} = \frac{1}{\lambda_{\min}(A_n)} = \frac{1}{4 \sin^2\left(\frac{\pi}{2(n+1)}\right)}

Poiché \sin(x) \approx x per x \to 0 e h = 1/(n+1):

\|A_n^{-1}\|_2 \approx \frac{1}{4 \left(\frac{\pi}{2(n+1)}\right)^2} = \frac{(n+1)^2}{\pi^2} = \frac{1}{\pi^2 h^2}

Quindi:

\|\mathbf{\epsilon}^{(n)}\|_2 \le h^2 \|A_n^{-1}\|_2 \|\mathbf{\tau}^{(n)}\|_\infty \approx h^2 \frac{1}{\pi^2 h^2} \|\mathbf{\tau}^{(n)}\|_\infty = \frac{1}{\pi^2} \|\mathbf{\tau}^{(n)}\|_\infty

Se la soluzione è regolare, \tau è limitato e la convergenza è garantita.

Analisi errore in norma infinito

Per stimare \|\mathbf{\epsilon}^{(n)}\|_\infty ci basta studiare \|A_n^{-1}\|_\infty. Notiamo che -A_n = \frac{2}{h^2}(I - B_n), dove B_n è una matrice tridiagonale con zeri sulla diagonale e 1/2 sopra e sotto. Per Gershgorin, gli autovalori di \text{trid}(1,0,1) sono tra 0 e 2. B_n ha raggio spettrale \rho(B_n) < 1, quindi la serie \sum B_n^i converge. Vale -A_n^{-1} = \sum B_n^i \ge 0. Segue inoltre che -A_n^{-1} > 0.

Analisi matriciale

Poiché -A_n^{-1} > 0, si ha \|A_n^{-1}\|_\infty = \|(-A_n^{-1}) \mathbf{e}\|_\infty, dove \mathbf{e} = (1, \dots, 1)^T. Ci interessa dunque trovare \mathbf{w} = -A_n^{-1} \mathbf{e}, ovvero risolvere A_n \mathbf{w} = -\mathbf{e}. Poiché A_n rappresenta la derivata seconda discretizzata, possiamo intuire che \mathbf{w} sia legato alla discretizzazione di una funzione u(x) tale che u''(x) = 1. La soluzione è u(x) = \frac{1}{2}x(1-x). I valori di w_i dovrebbero essere proporzionali a u(x_i) = \frac{1}{2} x_i(1-x_i). Infatti, si può verificare che w_i = \frac{h^2}{2} i(n+1-i). Con questa scelta, si ha:

(A_n \mathbf{w})_i = -\frac{1}{h^2} (w_{i-1} - 2w_i + w_{i+1}) = -1 = (-\mathbf{e})_i

Quindi:

\|A_n^{-1}\|_\infty = \|\mathbf{w}\|_\infty = \max_i \frac{h^2}{2} i(n+1-i)

Il massimo si ha per i \approx (n+1)/2:

\|A_n^{-1}\|_\infty \approx \frac{h^2}{2} \frac{n+1}{2} \left(n+1 - \frac{n+1}{2}\right) = \frac{h^2}{8}(n+1)^2 = \frac{1}{8}

Dunque, l’errore globale in norma infinito è:

\|\mathbf{\epsilon}^{(n)}\|_\infty \le h^2 \|A_n^{-1}\|_\infty \|\mathbf{\tau}^{(n)}\|_\infty \approx \frac{h^2}{8} \|\mathbf{\tau}^{(n)}\|_\infty

Se u \in C^4([0,1]), allora \|\mathbf{\tau}^{(n)}\|_\infty \le \frac{1}{12} \max_{[0,1]} |u^{(4)}|, e quindi:

\|\mathbf{\epsilon}^{(n)}\|_\infty \le \frac{h^2}{96} \max_{[0,1]} |u^{(4)}|

Questo mostra che il metodo converge con ordine 2.

Analisi tramite il principio del massimo

Un approccio alternativo per la stima dell’errore si basa sul principio del massimo.

L[u] = -\sum_{i,j=1}^d a_{ij}(x) \frac{\partial^2 u}{\partial x_i \partial x_j} - \sum_{i=1}^d b_i(x) \frac{\partial u}{\partial x_i}

\max_{x \in \overline{\Omega}} u(x) = \max_{x \in \partial \Omega} u(x)

\max_{i=1,\dots,n} v_i \le \max(v_0, v_{n+1})

-\frac{1}{h^2}(v_{k-1} - 2v_k + v_{k+1}) \ge 0 \implies 2v_k \le v_{k-1} + v_{k+1}

\|\mathbf{\epsilon}^{(n)}\|_\infty \le \frac{h^2}{8} \|\mathbf{\tau}^{(n)}\|_\infty

A_n \mathbf{z}^{\pm} = \pm A_n \mathbf{\epsilon}^{(n)} + h^2 \|\mathbf{\tau}^{(n)}\|_\infty A_n \mathbf{w} = \pm h^2 \mathbf{\tau}^{(n)} + h^2 \|\mathbf{\tau}^{(n)}\|_\infty \mathbf{e} \ge 0

\pm \epsilon_i + h^2 \|\mathbf{\tau}^{(n)}\|_\infty w_i \le 0 \implies |\epsilon_i| \le h^2 \|\mathbf{\tau}^{(n)}\|_\infty w_i

\|\mathbf{\epsilon}^{(n)}\|_\infty \le h^2 \|\mathbf{\tau}^{(n)}\|_\infty \|\mathbf{w}\|_\infty

Caso multidimensionale

Siano u: \Omega \subseteq \mathbb{R}^d \to \mathbb{R} e f: \Omega \to \mathbb{R}, g: \partial \Omega \to \mathbb{R}.

\begin{cases} \Delta u = f & \text{su } \Omega \\ u = g & \text{su } \partial \Omega \end{cases}

consideriamo d=2 ed una discretizzazione su una griglia m \times n:

x_i = i h_x, \quad h_x = \frac{1}{m+1}

y_j = j h_y, \quad h_y = \frac{1}{n+1}

Consideriamo le approssimazioni per le derivate parziali:

\begin{align}
\frac{\partial^2 u}{\partial x^2} \bigg|_{x_i, y_j} &= \frac{1}{h_x^2} (u_{i-1,j} - 2u_{i,j} + u_{i+1,j}) + h_x^2 \tau_{ij} \\[1em]
\frac{\partial^2 u}{\partial y^2} \bigg|_{x_i, y_j} &= \frac{1}{h_y^2} (u_{i,j-1} - 2 u_{i,j} + u_{i,j+1}) + h_y^2 \nu_{ij}
\end{align}

O come STENCIL:

\begin{bmatrix} & -1 & \\ -1 & 4 & -1 \\ & -1 & \end{bmatrix}

Una possibile scelta per L_{m,n}[u] è con:

• L_m: derivata lungo i
• L_n: derivata lungo j

\rightsquigarrow L_{m,n}[u] = L_m U + U L_n

con U = (u_{ij})_{ij} matrice che rappresenta i valori di u(x_i, y_j). Il sistema sarà:

A_m U + U A_n = F + h_x^2 \tau^{(m,n)} + h_y^2 \nu^{(m,n)}

\rightsquigarrow A_m V + V A_n = F \quad (V \approx U)

Questo non è molto comodo da risolvere, proviamo a riordinarlo:

\text{vec}(U) := \begin{bmatrix} \mathbf{u}_{:,1} \\ \mathbf{u}_{:,2} \\ \vdots \\ \mathbf{u}_{:,n} \end{bmatrix}

Introduciamo il prodotto di Kronecker:

A \otimes B = \begin{bmatrix} a_{11} B & \dots & a_{1m} B \\ \vdots & & \vdots \\ a_{n1} B & \dots & a_{nm} B \end{bmatrix}

per A \in \mathbb{R}^{n_1 \times m_1} e B \in \mathbb{R}^{n_2 \times m_2}, A \otimes B sarà (n_1 n_2) \times (m_1 m_2).

Possiamo anche definire per V \in \mathbb{R}^n ssp, W \in \mathbb{R}^m ssp, V \otimes W come ssp in \mathbb{R}^{nm}.

\text{vec}(AXB) = (B^T \otimes A) \text{vec}(X)

Possiamo riscrivere il sistema precedente come

A_m V I_n + I_m V A_n = F

e applicando \text{vec}(\cdot):

(I_n \otimes A_m + A_n \otimes I_m) \text{vec}(V) = \text{vec}(F)

(A_n^T = A_n poiché simmetrica).

I sistemi della forma:

• Lyapunov: AX + XA^T = F
• Sylvester: AX + XB = F

Nel nostro caso per n=m possiamo ricondurci ad una matrice diagonale, ad esempio prendendo: A = S D S^T con S trasformata in seni e D sampling di \frac{1}{h^2}(2 - 2 \cos \theta).

E l’analisi dell’errore deriva da quella per il caso 1-D.

Parabolico

Ad esempio possiamo considerare:

\begin{cases}
\gamma \dfrac{\partial u}{\partial t} - \Delta u = 0 \\
u = u_0 & \text{su } \{0\} \times \Omega \\
u = g & \text{su } [0,T] \times \partial \Omega
\end{cases}

Se ad esempio discretizziamo u come u_i^j con u_i^j \approx u(t_j, x_i):

\rightsquigarrow \gamma \frac{u_i^{j+1} - u_i^j}{\Delta t} - \frac{1}{h^2}(u_{i-1}^j - 2 u_i^j + u_{i+1}^j) = h^2 \tau_{i}^j + \Delta t \dots

Da cui otteniamo la seguente senza termine di errore:

\gamma \frac{u_i^{j+1} - u_i^j}{\Delta t} - \frac{1}{h^2}(u_{i-1}^j - 2 u_i^j + u_{i+1}^j) = 0

\implies \mathbf{v}^{(j+1)} = \mathbf{v}^{(j)} + \frac{\Delta t}{\gamma h^2} A \mathbf{v}^{(j)} = \left( I + \frac{\Delta t}{\gamma h^2} A \right) \mathbf{v}^{(j)}

Possiamo scrivere il sistema per tutti i passi temporali come:

A_{m,n} \begin{bmatrix} \mathbf{v}^{(1)} \\ \vdots \\ \mathbf{v}^{(n)} \end{bmatrix} = \begin{bmatrix} \mathbf{v}^{(0)} \\ \vdots \end{bmatrix}

con A_{m,n} = \begin{bmatrix} I & & & \\ -T & I & & \\ & \ddots & \ddots & \\ & & -T & I \end{bmatrix} e T \coloneqq I + \frac{\Delta t}{\gamma h^2} A.

Bisogna stare attenti perché l’inversa di una matrice di questo tipo ha la forma:

\begin{bmatrix} I & & & \\ -T & I & & \\ & \ddots & \ddots & \\ & & -T & I \end{bmatrix}^{-1} = \begin{bmatrix} I & & & \\ T & I & & \\ T^2 & T & I & \\ \vdots & & & \ddots \end{bmatrix}

quindi bisogna avere \rho(T) < 1 per avere convergenza.