Integrazione numerica dell' equazione delle onde (1D)

Molti fenomeni della Fisica riconducibili a campi (grandezze continue definite in regioni dello spazio, eventualmente dipendenti dal tempo) sono descritti a partire da equazioni differenziali a derivate parziali (PDE). Equazioni cioè in cui sono date delle relazioni tra derivate parziali della/e funzione/i che descrivono il sistema, funzioni stesse ed eventualmente funzioni note delle variabili indipendenti.

Esempi di PDE sono:

• Equazione di Poisson (potenziale gravitazionale newtowniano, potenziale elettrostatico):
  
  

  $$\frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \p... ...ial y^2} + \frac{\partial^2 \phi}{\partial z^2} = \rho(x,y,z)$$ (1)

  
  

• Equazione di diffusione (variazioni spazio-temporali della temperatura, concentrazione):
  
  

  $$\frac{\partial u}{\partial t} = h^2 \left( \frac{\partial^2... ... u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2} \right)$$ (2)

  
  

• Equazione delle onde (onde elastiche, sonore, etc. in mezzi non dispersivi):
  
  

  $$v^2 \left( \frac{\partial^2 u}{\partial x^2} + \frac{\part... ...partial z^2} \right) - \frac{\partial^2 u}{\partial t^2} = 0$$ (3)

  
  

Per ciascuna delle equazioni precedenti, l' analisi fisico-matematica del problema suggerisce quali siano le condizioni più appropriate da richiedere alla funzione incognita per avere un problema matematicamente ben posto.

Qui ci occuperemo della soluzione numerica dell' equazione (3) nel caso particolare di funzioni che dipendano, oltre che dal tempo, da un' unica variabile spaziale. L' equazione d' onda unidimensionale:

$$\frac{\partial^2 u}{\partial t^2} = v^2 \frac{\partial^2 u}{\partial x^2}$$ (4)

soggetta ad opportune condizioni al contorno (restrizioni legate al dominio di $u(x,t)$ come funzione di $x$) e iniziali.

In particolare considereremo i seguenti casi di condizioni al contorno

• estremi fissi:
  

  $$u(-L,t) = u(L,t) = 0$$ (5)

  
  

• condizioni periodiche
  

  $$u(-L,t) = u(L,t)$$ (6)

  $$\frac{\partial u}{\partial x}(-L,t) = \frac{\partial u}{\partial x}(L,t)$$ (7)

  
  

Le condizioni iniziali sono specificate dando il campo $u(x,0)$ e la sua derivata parziale rispetto al tempo $\frac{\partial u}{\partial t}$ al tempo $t=0$ su tutto il dominio (rispettando le condizioni al contorno):

$$u(x,t=0) = f(x)$$ (8)

$$\frac{\partial u}{\partial t}(x,t=0) = g(x)$$ (9)

Un modo per risolvere numericamente l' equazione (4) è di usare un' approssimazione discretizzata per la derivata seconda rispetto a x.

Per esempio, usando la formula alle differenze simmetriche per la derivata prima numerica abbiamo:

$$f'(x) = \frac{f(x+h/2)-f(x-h/2)}{h}$$ (10)

a meno di termini $O(h^2)$. Quindi per la derivata seconda varrà:

$$f''(x) = \frac{f'(x+h/2)-f'(x-h/2)}{h} = \frac{f(x+h)-2 f(x) + f(x-h)}{h^2}$$ (11)

sempre a meno di termini $O(h^2)$.

Se discretizziamo quindi le funzioni di $x$ su una griglia spaziale di punti equispaziati $x_i$, possiamo indicare con $f_i(t)$ il valore di una qualsiasi funzione $f(x,t)$ nell' i-esimo punto della griglia.

L' equazione (4) viene allora trasformata dalla discretizzazione spaziale in un sistema di tante equazioni differenzali ordinarie (nel tempo) quanti sono i punti indipendenti della griglia spaziale:

$$\frac{d^2 u_i}{d t^2} = \left(\frac{v}{h}\right)^2 (u_{i+1}(t) - 2 u_i(t) + u_{i-1}(t)).$$ (12)

Se gli estremi sono fissi e l' indice $i$ varia tra $-N$ e $N$, occorrerà imporre

$$u_i = 0$$ (13)

$$\dot u_i = 0$$ (14)

per $i=-N$ e $i=N$. Nel caso di condizioni periodiche avremo:

$$u_N = u_{-N}$$ (15)

$$\dot u_N = \dot u_{-N}.$$ (16)

Per la soluzione numerica del sistema (12) possiamo utilizzare l' algoritmo di Verlet, già utiizzato per la dinamica del punto.

A differenza del caso della dinamica del punto, nel caso dell' equazione d' onda abbiamo due parametri numerici, l' intervallino di discretizzazione spaziale ($h$) e quello di discretizzazione temporale ($\Delta t$) che vanno fissati e da cui dipenderà l' accuratezza dell' integrazione delle equazioni del moto.

Il vincolo su $h$ è che permetta di rappresentare in modo soddisfacente la funzione e le sue derivate spaziali.

Una volta scelto $h$, la discretizzazione della derivata seconda spaziale e l' algoritmo di integrazione numerica, il valore di $\Delta t$ non è completamente indipendente ma dovrà rispettare la cosiddetta condizione di Courant:

$$\Delta t \leq C \frac{h}{v}$$

dove $C$ è una costante (vale $1$ per l' algoritmo di Verlet).

Come al solito nel caso dell' integrazione numerica di equazioni del moto conservative, la verifica della qualità della conservazione dell' energia darà una misura della qualità dell' integrazione.