# Metodo di Laplace

Nell'analisi matematica, il metodo di Laplace, il cui nome deriva da Pierre-Simon Laplace, è una tecnica usata per approssimare integrali nella forma
    $\int _{a}^{b}\!e^{Mf(x)}\,dx,$
dove $f(x)$ è una qualunque funzione derivabile due volte, $M$ è un numero "grande" e gli estremi d'integrazione $a$ e $b$ possono essere anche infiniti. Questa tecnica fu per la prima volta presentata nell'articolo "Mémoir sur la probabilité des causes par évènemens" di Laplace del 1774.

## L'idea del metodo di Laplace

Si assuma che la funzione $f(x)$ abbia un unico massimo globale in $x_{0}$. Allora, il valore $f(x_{0})$ sarà più grande degli altri valori di $f(x)$. Se si moltiplica questa funzione per un numero grande $M$, il rapporto fra $Mf(x_{0})$ e $Mf(x)$ rimane lo stesso (poiché $Mf(x_{0})/Mf(x)=f(x_{0})/f(x))$, ma crescerà esponenzialmente nella funzione $e^{Mf(x)}$ (vedere figura). Perciò solo punti $x$ in un intorno di $x_{0}$ daranno significativi contributi all'integrale della funzione, che può essere stimato.
Per descrivere e motivare il metodo, sono necessarie alcune ipotesi. Si assuma che $x_{0}$ non sia un estremo di integrazione, che il valore di $f(x)$ non possa essere molto vicino a $f(x_{0})$ a meno che $x$ sia vicino a $x_{0}$ e che la derivata seconda $f''(x_{0})<0$.
Si può sviluppare $f(x)$ intorno a $x_{0}$ usando il teorema di Taylor e ottenendo
    $f(x)=f(x_{0})+f'(x_{0})(x-x_{0})+{\frac {1}{2}}f''(x_{0})(x-x_{0})^{2}+R,$
dove $R=O\left((x-x_{0})^{3}\right).$
Poiché $f$ ha un massimo globale in $x_{0}$, e siccome $x_{0}$ non è un estremo, esso è un punto stazionario, perciò la derivata di $f$ in $x_{0}$ si annulla. Dunque, la funzione può essere approssimata al secondo ordine come
    $f(x)\approx f(x_{0})-{\frac {1}{2}}|f''(x_{0})|(x-x_{0})^{2},$
per $x$ vicino a $x_{0}$ (si ricordi che la derivata seconda nel punto di massimo $x_{0}$ è negativa). Le ipotesi assicurano la precisione dell'approssimazione
    $\int _{a}^{b}\!e^{Mf(x)}\,dx\approx e^{Mf(x_{0})}\int _{a}^{b}e^{-M|f''(x_{0})|(x-x_{0})^{2}/2}\,dx$
(vedere la figura sulla destra). Quest'ultimo integrale sarebbe un integrale di Gauss se i limiti di integrazione andassero da $-\infty$ a $+\infty$ (che può essere assunto dato che l'esponenziale decade molto velocemente lontano da $x_{0}$), e pertanto può essere calcolato. Si trova così che
    $\int _{a}^{b}\!e^{Mf(x)}\,dx\approx {\sqrt {\frac {2\pi }{M|f''(x_{0})|}}}e^{Mf(x_{0})},\qquad {\text{con }}M\to +\infty .$

## Enunciato

Si assuma che $f(x)$ sia una funzione di classe $C^{2}$ su $[a,b]$ con $x_{0}\in (a,b)$ l'unico punto tale che $f(x_{0})=\max _{[a,b]}f(x)$. Si assuma inoltre che $f''(x_{0})<0$.
Allora
    $\lim _{n\to +\infty }{\frac {\int _{a}^{b}e^{nf(x)}\,dx}{e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}}}=1.$

## Dimostrazione

Minorante:
Sia $\varepsilon >0$. Per la continuità di $f''$ esiste $\delta >0$ tale che se $|x_{0}-c|<\delta$ allora $f''(c)\geq f''(x_{0})-\varepsilon$. Dal teorema di Taylor, per ogni $x\in (x_{0}-\delta ,x_{0}+\delta )$, $f(x)\geq f(x_{0})+{\frac {1}{2}}(f''(x_{0})-\varepsilon )(x-x_{0})^{2}$.
Quindi si ha il seguente minorante:
    $\int _{a}^{b}e^{nf(x)}\,dx\geq \int _{x_{0}-\delta }^{x_{0}+\delta }e^{nf(x)}\,dx\geq e^{nf(x_{0})}\int _{x_{0}-\delta }^{x_{0}+\delta }e^{{\frac {n}{2}}(f''(x_{0})-\varepsilon )(x-x_{0})^{2}}\,dx=e^{nf(x_{0})}{\sqrt {\frac {1}{n(-f''(x_{0})+\varepsilon )}}}\cdot \int _{-\delta {\sqrt {n(-f''(x_{0})+\varepsilon )}}}^{\delta {\sqrt {n(-f''(x_{0})+\varepsilon )}}}e^{-{\frac {1}{2}}y^{2}}\,dy,$
dove l'ultima uguaglianza è stata ottenuta dal cambio di variabili $y={\sqrt {n(-f''(x_{0})+\varepsilon )}}(x-x_{0})$. Si ricordi che $f''(x_{0})<0$ e quindi è possibile estrarne la radice quadrata.
Se si dividono entrambi i membri della precedente disuguaglianza per $e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}$ e se ne prende il limite si ottiene:
    $\lim _{n\to +\infty }{\frac {\int _{a}^{b}e^{nf(x)}\,dx}{e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}}}\geq \lim _{n\to +\infty }{\frac {1}{\sqrt {2\pi }}}\int _{-\delta {\sqrt {n(-f''(x_{0})+\varepsilon )}}}^{\delta {\sqrt {n(-f''(x_{0})+\varepsilon )}}}e^{-{\frac {1}{2}}y^{2}}\,dy\,\cdot {\sqrt {\frac {-f''(x_{0})}{-f''(x_{0})+\varepsilon }}}={\sqrt {\frac {-f''(x_{0})}{-f''(x_{0})+\varepsilon }}}.$
Poiché è vero per un arbitrario $\varepsilon$, si trova il minorante:
    $\lim _{n\to +\infty }{\frac {\int _{a}^{b}e^{nf(x)}\,dx}{e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}}}\geq 1.$
Da notare che la dimostrazione funziona anche quando $a=-\infty$ oppure $b=+\infty$ (o entrambi).
Maggiorante: La dimostrazione del maggiorante è simile a quella del minorante ma ci sono alcuni inconvenienti. Di nuovo si inizia prendendo un $\varepsilon >0$ ma allo scopo di far funzionare la dimostrazione serve che $\varepsilon$ sia piccolo abbastanza affinché $f''(x_{0})+\varepsilon <0$. Quindi, come sopra, dalla continuità di $f''$ e il teorema di Taylor si trova $\delta >0$ tale che se $|x-x_{0}|<\delta$, allora $f(x)\leq f(x_{0})+{\frac {1}{2}}(f''(x_{0})+\varepsilon )(x-x_{0})^{2}$. Infine per ipotesi (assumendo $a,b$ finiti) esiste un $\eta >0$ tale che se $|x-x_{0}|\geq \delta$, allora $f(x)\leq f(x_{0})-\eta$.
Si può ora calcolare il seguente maggiorante:
    ${\begin{aligned}&\int _{a}^{b}e^{nf(x)}\,dx\leq \int _{a}^{x_{0}-\delta }e^{nf(x)}\,dx+\int _{x_{0}-\delta }^{x_{0}+\delta }e^{nf(x)}\,dx+\int _{x_{0}+\delta }^{b}e^{nf(x)}\,dx\leq (b-a)e^{n(f(x_{0})-\eta )}+\int _{x_{0}-\delta }^{x_{0}+\delta }e^{nf(x)}\,dx\\&\leq (b-a)e^{n(f(x_{0})-\eta )}+e^{nf(x_{0})}\int _{x_{0}-\delta }^{x_{0}+\delta }e^{{\frac {n}{2}}(f''(x_{0})+\varepsilon )(x-x_{0})^{2}}\,dx\leq (b-a)e^{n(f(x_{0})-\eta )}+e^{nf(x_{0})}\int _{-\infty }^{+\infty }e^{{\frac {n}{2}}(f''(x_{0})+\varepsilon )(x-x_{0})^{2}}\,dx\\&\leq (b-a)e^{n(f(x_{0})-\eta )}+e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0})-\varepsilon )}}}.\end{aligned}}$
Se si dividono entrambi i membri della disuguaglianza per $e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}$ e se ne prende il limite si ottiene:
    $\lim _{n\to +\infty }{\frac {\int _{a}^{b}e^{nf(x)}\,dx}{e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}}}\leq \lim _{n\to +\infty }(b-a)e^{-\eta n}{\sqrt {\frac {n(-f''(x_{0}))}{2\pi }}}+{\sqrt {\frac {-f''(x_{0})}{-f''(x_{0})-\varepsilon }}}={\sqrt {\frac {-f''(x_{0})}{-f''(x_{0})-\varepsilon }}}.$
Poiché $\varepsilon$ è arbitrario si ha il maggiorante:
    $\lim _{n\to +\infty }{\frac {\int _{a}^{b}e^{nf(x)}\,dx}{e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}}}\leq 1.$
E combinandolo con il risultato ricavato prima si dimostra l'enunciato.
Da notare che la dimostrazione precedente fallisce quando $a=-\infty$ oppure $b=+\infty$ (o entrambi). Per trattare questi casi, c'è bisogno di ulteriori ipotesi. Un'assunzione sufficiente (e non necessaria) è che per $n=1$, l'integrale $\int _{a}^{b}e^{nf(x)}\,dx$ sia finito, e che il numero $\eta$ come sopra esista (si osservi che questa deve essere un'ipotesi solo nel caso di $a$ o $b$ infinito). La dimostrazione procede altrimenti come prima, ma gli integrali
    $\int _{a}^{x_{0}-\delta }e^{nf(x)}\,dx+\int _{x_{0}+\delta }^{b}e^{nf(x)}\,dx$
devono essere stimati superiormente da
    $\int _{a}^{x_{0}-\delta }e^{nf(x)}\,dx+\int _{x_{0}+\delta }^{b}e^{nf(x)}\,dx\leq \int _{a}^{b}e^{f(x)}e^{(n-1)(f(x_{0})-\eta )}\,dx=e^{(n-1)(f(x_{0})-\eta )}\int _{a}^{b}e^{f(x)}\,dx$
invece di $(b-a)e^{n(f(x_{0})-\eta )}$ come per il minorante, così che quando si divide per $e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}$, si ottiene per questo termine
    ${\frac {e^{(n-1)(f(x_{0})-\eta )}\int _{a}^{b}e^{f(x)}\,dx}{e^{nf(x_{0})}{\sqrt {\frac {2\pi }{n(-f''(x_{0}))}}}}}=e^{-(n-1)\eta }{\sqrt {n}}e^{-f(x_{0})}\int _{a}^{b}e^{f(x)}\,dx{\sqrt {\frac {-f''(x_{0})}{2\pi }}}$
il cui limite per $n\to +\infty$ è $0$. Il resto della dimostrazione (l'analisi dei termini dominanti) procede come sopra.
La condizione data nel caso di intervallo infinito è, come detto precedentemente, sufficiente ma non necessaria. Tuttavia, la condizione è soddisfatta nella maggior parte delle applicazioni: la condizione semplicemente afferma che l'integrale che si sta studiando sia ben definito (non infinito) e che il massimo della funzione in $x_{0}$ sia un "vero" massimo (il numero $\eta >0$ deve esistere). Non c'è inoltre bisogno di richiedere che l'integrale sia finito per $n=1$ ma è sufficiente che lo sia per un qualche $n=N$.

## Applicazione: approssimazione di Stirling

Il metodo di Laplace può essere utilizzato per derivare l'approssimazione di Stirling
    $N!\approx {\sqrt {2\pi N}}N^{N}e^{-N},$
per un intero $N$ grande.
Dalla definizione della funzione Gamma, si ha
    $N!=\Gamma (N+1)=\int _{0}^{\infty }e^{-x}x^{N}\,dx.$
Ora effettuando il cambio di variabile
    $x=Nz$
si ottiene
    $N!=\int _{0}^{\infty }e^{-Nz}\left(Nz\right)^{N}N\,dz=N^{N+1}\int _{0}^{\infty }e^{-Nz}z^{N}\,dz=N^{N+1}\int _{0}^{\infty }e^{-Nz}e^{N\ln z}\,dz=N^{N+1}\int _{0}^{\infty }e^{N(\ln z-z)}\,dz.$
Questo integrale ha la forma necessaria per il metodo di Laplace con
    $f\left(z\right)=\ln {z}-z$
che è derivabile con continuità due volte:
    $f'(z)={\frac {1}{z}}-1,$
    $f''(z)=-{\frac {1}{z^{2}}}.$
Il massimo di $f(x)$ si trova $z_{0}=1$, e la derivata seconda in quel punto ha valore $-1$. Pertanto, si ricava
    $N!\approx N^{N+1}{\sqrt {\frac {2\pi }{N}}}e^{-N}={\sqrt {2\pi N}}N^{N}e^{-N}.$

## Generalizzazioni

L'approssimazione di Laplace può essere generalizzata agli integrali nella forma
    $\int _{a}^{b}h(x)e^{Mg(x)}\,dx\approx {\sqrt {\frac {2\pi }{M|g''(x_{0})|}}}h(x_{0})e^{Mg(x_{0})},\qquad {\text{con }}M\to +\infty ,$
dove $h$ è positiva. È importante sottolineare che la precisione dell'approssimazione dipende dalla variabile di integrazione.
Nel caso a più variabili, dove $\mathbf {x}$ è un vettore $n$-dimensionale e $f(\mathbf {x} )$ è una funzione scalare di $\mathbf {x}$, l'approssimazione di Laplace diventa:
    $\int e^{Mf(\mathbf {x} )}\,d\mathbf {x} \approx \left({\frac {2\pi }{M}}\right)^{n/2}|-H(f)(\mathbf {x} _{0})|^{-1/2}e^{Mf(\mathbf {x} _{0})},\qquad {\text{con }}M\to +\infty ,$
con $H(f)(\mathbf {x} _{0})$ la matrice hessiana di $f$ calcolata in $\mathbf {x} _{0}$ e dove $|\cdot |$ indica il determinante. Analogamente al caso di una variabile, la matrice hessiana deve essere definita negativa.

## Precisione del metodo

Prima di tutto, si ponga senza perdita di generalità che il massimo globale si trovi in $x_{0}=0$. Perciò, quello che si vuole è l'errore relativo $\left|R\right|$ come mostrato sotto
    $\int _{a}^{b}\!h(x)e^{Mg(x)}\,dx=h(0)e^{Mg(0)}s\underbrace {\int _{a/s}^{b/s}{\frac {h(x)}{h(0)}}e^{M\left[g(sy)-g(0)\right]}dy} _{1+R},$
dove $s\equiv {\sqrt {\frac {2\pi }{M\left|g''(0)\right|}}}$. Quindi, posto $A\equiv {\frac {h(sy)}{h(0)}}e^{M\left[g(sy)-g(0)\right]}$ e $A_{0}\equiv e^{-\pi y^{2}}$, si ottiene
    $\left|R\right|=\left|\int _{a/s}^{b/s}A\,dy-\int _{-\infty }^{\infty }A_{0}\,dy\right|$
poiché $\int _{-\infty }^{\infty }A_{0}\,dy=1$. Ora si deve trovare un maggiorante.
Grazie a $\left|A+B\right|\leq |A|+|B|$, si può separare l'integrazione in 5 parti di 3 differenti tipi: $(a)$, $(b)$ e $(c)$, rispettivamente. Pertanto,
    $|R|<\underbrace {\left|\int _{D_{y}}^{\infty }A_{0}dy\right|} _{(a_{1})}+\underbrace {\left|\int _{D_{y}}^{b/s}Ady\right|} _{(b_{1})}+\underbrace {\left|\int _{-D_{y}}^{D_{y}}\left(A-A_{0}\right)dy\right|} _{(c)}+\underbrace {\left|\int _{a/s}^{-D_{y}}Ady\right|} _{(b_{2})}+\underbrace {\left|\int _{-\infty }^{-D_{y}}A_{0}dy\right|} _{(a_{2})},$
dove $(a_{1})$ e $(a_{2})$ sono simili, quindi si calcolerà solo $(a_{1})$, e analogamente per $(b_{1})$ e $(b_{2})$.
Per $(a_{1})$, dopo aver rinominato $z\equiv \pi y^{2}$, si ha
    $(a_{1})=\left|{\frac {1}{2{\sqrt {\pi }}}}\int _{\pi D_{y}^{2}}^{\infty }e^{-z}z^{-1/2}dz\right|<{\frac {e^{-\pi D_{y}^{2}}}{2\pi D_{y}}}.$
Questo significa che fintanto che $D_{y}$ è abbastanza grande, esso tenderà a zero.
Per $(b_{1})$, si ricava
    $(b_{1})\leq \left|\int _{D_{y}}^{b/s}\left[{\frac {h(sy)}{h(0)}}\right]_{\mathrm {max} }e^{Mm(sy)}dy\right|,$
dove
    $m(x)\geq g(x)-g(0),\qquad {\text{con }}x\in [sD_{y},b]$
e $h(x)$ dovrebbe avere lo stesso segno di $h(0)$ nella regione. Si scelga $m(x)$ come la tangente in $x=sD_{y}$, cioè $m(sy)=g(sD_{y})-g(0)+g'(sD_{y})\left(sy-sD_{y}\right)$ (che è mostrata in figura).
Dalla figura si può osservare che quando $s$ o $D_{y}$ diventa piccolo, la regione che soddisfa la precedente disuguaglianza diventa sempre più grande. Dunque, se si vuole trovare una $m(x)$ adatta a coprire l'intera $f(x)$ nell'intervallo di $(b_{1})$, $D_{y}$ deve avere un estremo superiore. Inoltre, siccome l'integrale $e^{-\alpha x}$ è semplice, si userà per stimare l'errore relativo dovuto a $(b_{1})$.
Usando lo sviluppo di Taylor, si ottiene
    $M\left[g(sD_{y})-g(0)\right]=M\left[{\frac {g''(0)}{2}}s^{2}D_{y}^{2}+{\frac {g'''(\xi )}{6}}s^{3}D_{y}^{3}\right]\,\,{\text{con}}\,\,\xi \in [0,sD_{y}]=-\pi D_{y}^{2}+{\frac {(2\pi )^{3/2}g'''(\xi )D_{y}^{3}}{6{\sqrt {M}}|g''(0)|^{3/2}}},$
e
    $Msg'(sD_{y})=Ms\left(g''(0)sD_{y}+{\frac {g'''(\zeta )}{2}}s^{2}D_{y}^{2}\right),\qquad {\text{con }}\zeta \in [0,sD_{y}]=-2\pi D_{y}+{\sqrt {\frac {2}{M}}}\left({\frac {\pi }{|g''(0)|}}\right)^{3/2}g'''(\zeta )D_{y}^{2},$
e si sostituisce nel calcolo di $(b_{1})$. Tuttavia, si trova che i resti dei due sviluppi sono entrambi inversamente proporzionali alla radice di $M$, perciò si tralasceranno per rendere più elegante il calcolo.
    $(b_{1})\leq \left|\left[{\frac {h(sy)}{h(0)}}\right]_{\text{max}}e^{-\pi D_{y}^{2}}\int _{0}^{b/s-D_{y}}e^{-2\pi D_{y}y}dy\right|\leq \left|\left[{\frac {h(sy)}{h(0)}}\right]_{\text{max}}e^{-\pi D_{y}^{2}}{\frac {1}{2\pi D_{y}}}\right|.$
In aggiunta, tenderà a zero quando $D_{y}$ diventa arbitrariamente grande, ma non si dimentichi che il limite superiore di $D_{y}$ deve essere considerato nel calcolo.
A proposito dell'integrazione vicino a $x=0$, si può usare anche il teorema di Taylor per calcolarlo. Quando $h'(0)\neq 0$
    $(c)\leq \int _{-D_{y}}^{D_{y}}e^{-\pi y^{2}}\left|{\frac {sh'(\xi )}{h(0)}}y\right|\,dy<{\sqrt {\frac {2}{\pi M|g''(0)|}}}\left|{\frac {h'(\xi )}{h(0)}}\right|_{\text{max}}\left(1-e^{-\pi D_{y}^{2}}\right)$
e si trova che è inversamente proporzionale a ${\sqrt {M}}$. Infatti, $(c)$ avrà lo stesso comportamento quando $h(x)$ è costante.
Infine, l'integrazione vicino al punto stazionario diventa piccola quando ${\sqrt {M}}$ diventa grande, e le parti rimanenti tenderanno a zero fintanto che $D_{y}$ è abbastanza grande, ma quest'ultimo ha un estremo superiore dovuto alla condizione che la funzione $m(x)$ è sempre maggiore di $g(x)-g(0)$ nella regione rimanente. Tuttavia, fino a che si trova un $m(x)$ che soddisfa la condizione, l'estremo superiore di $D_{y}$ può essere scelto direttamente proporzionale a ${\sqrt {M}}$ poiché $m(x)$ è la tangente di $g(x)-g(0)$ in $x=sD_{y}$. Quindi, più grande è $M$, più può essere grande $D_{y}$.

## Estensione del metodo di Laplace: la discesa del gradiente

Un'estensione del metodo di Laplace all'analisi complessa, insieme alla formula integrale di Cauchy, è usata per trovare un contorno di "discesa più ripida" per un (asintoticamente per grandi $M$) integrale equivalente, espresso come un integrale di linea. In particolare, se non esistono punti sulla retta reale in cui la derivata di $f$ si annulla, può essere necessario deformare in contorno di integrazione in uno ottimale, dove l'analisi discussa prima è possibile. Ancora l'idea principale è di ridurre, almeno in modo asintotico, il calcolo del dato integrale a uno più semplice e che quindi può essere valutato esplicitamente. Si veda il libro di Erdelyi (1956) per una semplice discussione (dove il metodo è chiamato "discesa del gradiente")
L'appropriata formulazione per il piano complesso è
    $\int _{a}^{b}\!e^{Mf(z)}\,dz\approx {\sqrt {\frac {2\pi }{-Mf''(z_{0})}}}e^{Mf(z_{0})},\qquad {\text{con }}M\to +\infty .$
per un percorso passante attraverso il punto di sella in $z_{0}$. Da notare l'esplicita presenza di un segno meno ad indicare la direzione della derivata seconda: non se ne può prendere il modulo. Inoltre se la funzione integranda è meromorfa, si può dover aggiungere i residui corrispondenti ai poli attraversati durante la deformazione del contorno (vedere per esempio la sezione 3 dell'artico di Okounkov "Symmetric functions and random partitions").