Metodi numerici per ODE basati sulla serie di Taylor
Appunti di calcolo numerico
Consideriamo la seguente equazione differenziale ordinaria (ODE):
\left\{ \begin{aligned} y' &= f(x,y) \\ y(a) &= y_a \end{aligned} \right. \tag{1}
Sotto le usuali condizioni di regolarità per f(x, y), possiamo garantire l’esistenza e l’unicità della soluzione.
Per approssimare numericamente la soluzione nell’intervallo [a, b], dividiamo l’intervallo in n parti uguali con ampiezza h = (b - a) / n. Cerchiamo di calcolare la soluzione in corrispondenza dei punti:
x_k = a + kh, \qquad k=0,1,\dots,n.
La soluzione esatta y(x) soddisfa:
y'(x) = f(x,y(x)), \qquad y(a) = y_a.
Possiamo sfruttare lo sviluppo in serie di Taylor per costruire metodi numerici che approssimano la soluzione.
Metodo di Eulero e Serie di Taylor
Utilizziamo la serie di Taylor per approssimare la soluzione y(x) al primo ordine:
y(x_{k+1}) = y(x_k) + y'(x_k) h + \frac{y''(\zeta)}{2}h^2, \qquad \zeta \in (x_k, x_{k+1}),
dove y'(x_k) si ottiene direttamente dalla ODE come y'(x_k) = f(x_k, y(x_k)). Otteniamo così:
y(x_{k+1}) = y(x_k) + f(x_k, y(x_k)) h + \frac{y''(\zeta)}{2}h^2.
Trascurando il termine d’errore \mathcal{O}(h^2), otteniamo il metodo di Eulero esplicito:
y_{k+1} = y_k + f(x_k, y_k) h,
che è un metodo del primo ordine.
Metodo di Taylor al Secondo Ordine
Per ottenere una migliore approssimazione, possiamo considerare un’espansione della serie di Taylor fino al secondo ordine:
y(x_{k+1}) = y(x_k) + y'(x_k) h + \frac{y''(x_k)}{2} h^2 + \frac{y'''(\zeta)}{6} h^3,\quad \zeta\in(x_k,x_{k+1})
Il problema che emerge è il calcolo di y''(x_k), che non è immediatamente disponibile. Tuttavia, possiamo calcolarlo derivando la ODE rispetto a x:
y''(x) = \frac{\mathrm{d}}{\mathrm{d}x} y'(x) = \frac{\mathrm{d}}{\mathrm{d}x} f(x, y(x)).
Applicando la chain rule, otteniamo:
y''(x) = \frac{\partial f(x, y(x))}{\partial x} + \frac{\partial f(x, y(x))}{\partial y} \cdot y'(x).
Sostituendo y'(x) con f(x, y(x)), possiamo esprimere y''(x) come:
y''(x) = \frac{\partial f(x, y(x))}{\partial x} + \frac{\partial f(x, y(x))}{\partial y} \cdot f(x, y(x)).
Sostituendo questa espressione nello sviluppo di Taylor, otteniamo:
\begin{aligned} y(x_{k+1}) &= y(x_k) + f(x_k, y(x_k)) h \\ &\quad + \frac{h^2}{2} \left( \frac{\partial f(x_k, y(x_k))}{\partial x} + \frac{\partial f(x_k, y(x_k))}{\partial y} \cdot f(x_k, y(x_k)) \right) \\ &\quad + \mathcal{O}(h^3). \end{aligned}
Trascurando il termine \mathcal{O}(h^3), otteniamo un metodo di Taylor al secondo ordine:
y_{k+1} = y_k + f(x_k, y_k) h + \frac{h^2}{2} \left( \frac{\partial f(x_k, y_k)}{\partial x} + \frac{\partial f(x_k, y_k)}{\partial y} \cdot f(x_k, y_k) \right).
Osservazione 1. Il calcolo delle derivate parziali di f(x, y) è essenziale per costruire metodi di Taylor di ordine superiore. Tuttavia, per funzioni complesse, il calcolo manuale delle derivate può risultare molto laborioso e soggetto a errori. In questi casi, l’uso del calcolo simbolico diventa fondamentale.
Software come Mathematica, Maple, o librerie di Python come SymPy permettono di automatizzare il processo di differenziazione simbolica.
Metodo di Taylor al Terzo Ordine
Per ottenere una migliore approssimazione della soluzione, possiamo estendere la serie di Taylor fino al terzo ordine:
\begin{aligned} y(x_{k+1}) &= y(x_k) + y'(x_k) h + \frac{y''(x_k)}{2} h^2 \\ &\quad + \frac{y'''(x_k)}{6} h^3 + \frac{y^{(4)}(\zeta)}{24} h^4, \quad \zeta \in (x_k, x_{k+1}). \end{aligned}
Il calcolo di y''(x_k) è già stato fatto nel metodo di Taylor al secondo ordine. Per y'''(x_k), possiamo procedere come segue:
\begin{aligned} y'''(x) &= \frac{\mathrm{d}}{\mathrm{d}x} y''(x) \\ &= \frac{\mathrm{d}}{\mathrm{d}x} \left( \frac{\partial f(x, y(x))}{\partial x} + \frac{\partial f(x, y(x))}{\partial y} \cdot f(x, y(x)) \right) \\ &= \frac{\partial^2 f(x, y(x))}{\partial x^2} + \frac{\partial^2 f(x, y(x))}{\partial y^2} f(x, y(x))^2 \\ &\quad + \left( \frac{\partial f(x, y(x))}{\partial y} \right)^2 f(x, y(x)) + 2 \frac{\partial^2 f(x, y(x))}{\partial x \partial y} f(x, y(x)) \\ &\quad + \frac{\partial f(x, y(x))}{\partial x} \cdot \frac{\partial f(x, y(x))}{\partial y}. \end{aligned}
Notazione Compatta
Per semplificare la scrittura e il calcolo, possiamo introdurre alcune scorciatoie:
\begin{aligned} f_k &= f(x_k, y(x_k)), \\ f_{x,k} &= \frac{\partial f(x_k, y(x_k))}{\partial x}, \\ f_{y,k} &= \frac{\partial f(x_k, y(x_k))}{\partial y}, \\ f_{xx,k} &= \frac{\partial^2 f(x_k, y(x_k))}{\partial x^2}, \\ f_{yy,k} &= \frac{\partial^2 f(x_k, y(x_k))}{\partial y^2}, \\ f_{xy,k} &= \frac{\partial^2 f(x_k, y(x_k))}{\partial x \partial y}. \end{aligned}
Sostituendo queste espressioni nello sviluppo di Taylor, otteniamo:
\begin{aligned} y(x_{k+1}) &= y(x_k) + h f_k \\ &\quad + \frac{h^2}{2} \left( f_{x,k} + f_{y,k} f_k \right) \\ &\quad + \frac{h^3}{6} \left( f_{xx,k} + f_{yy,k} f_k^2 + ( f_{y,k}^2 + 2 f_{xy,k} ) f_k + f_{x,k} f_{y,k} \right) \\ &\quad + \mathcal{O}(h^4). \end{aligned}
Metodo di Taylor al Terzo Ordine
Trascurando il termine d’errore \mathcal{O}(h^4), otteniamo il metodo di Taylor al terzo ordine:
\begin{aligned} y_{k+1} &= y_k + f(x_k, y_k) h \\ &\quad + \frac{h^2}{2} \left( \frac{\partial f(x_k, y_k)}{\partial x} + \frac{\partial f(x_k, y_k)}{\partial y} \cdot f(x_k, y_k) \right) \\ &\quad + \frac{h^3}{6} \left( f_{xx,k} + f_{yy,k} f_k^2 + ( f_{y,k}^2 + 2 f_{xy,k} ) f_k + f_{x,k} f_{y,k} \right). \end{aligned}
Considerazioni sul metodo
I metodi basati su Taylor richiedono il calcolo di derivate parziali multiple della funzione f(x, y). Questo pone alcune sfide significative:
Derivate simboliche: Per usare un metodo di ordine superiore, è necessario calcolare simbolicamente (o con un software di algebra computazionale) un numero crescente di derivate parziali.
Crescita della complessità: La complessità del metodo cresce rapidamente con l’ordine, rendendo l’implementazione e il calcolo più difficili per ordini superiori al terzo.
Tuttavia, se la forma analitica di f(x, y) è nota in anticipo, è possibile costruire un metodo numerico basato sulla serie di Taylor specifico per quella funzione, ottenendo così un’accurata approssimazione della soluzione.
Metodo basato su Taylor con f(x,y) nota
Sia y(x) la soluzione analitica della seguente equazione differenziale ordinaria (ODE):
y'(x) = f(x, y(x)).
Consideriamo la serie di Taylor della soluzione:
\begin{aligned} y(x+h) &= y(x) + h y'(x) + \frac{h^2}{2} y''(x) + \cdots + \frac{h^p}{p!} y^{(p)}(x) + \mathcal{O}(h^{p+1}), \end{aligned} \tag{2}
dove ogni derivata k-esima può essere espressa in termini della funzione f(x, y(x)). A partire da questa relazione, possiamo costruire delle funzioni D_k(x, y) che hanno la proprietà:
y^{(k)}(x) = D_k(x, y(x)),
ovvero, se valutate sulla soluzione esatta y(x), restituiscono la derivata k-esima della soluzione stessa.
Definizione Ricorsiva delle Funzioni D_k(x, y)
La costruzione delle funzioni D_k(x, y) è ricorsiva e parte dalla definizione:
D_1(x, y) = f(x, y),
poiché y'(x) = f(x, y(x)) = D_1(x, y(x)).
Se supponiamo che y^{(k)}(x) = D_k(x, y(x)), possiamo derivare y^{(k+1)}(x) come segue:
\begin{aligned} y^{(k+1)}(x) &= \frac{\mathrm{d}}{\mathrm{d}x} y^{(k)}(x) \\ &= \frac{\mathrm{d}}{\mathrm{d}x} D_k(x, y(x)) \\ &= \frac{\partial D_k(x, y(x))}{\partial x} + \frac{\partial D_k(x, y(x))}{\partial y} \cdot y'(x) \\ &= \frac{\partial D_k(x, y(x))}{\partial x} + \frac{\partial D_k(x, y(x))}{\partial y} \cdot f(x, y(x)). \end{aligned}
Da questa relazione, otteniamo la seguente formula ricorsiva per D_{k+1}(x, y):
\left[ \begin{aligned} D_1(x, y) &= f(x, y), \\ D_{k+1}(x, y) &= \frac{\partial D_k(x, y)}{\partial x} + \frac{\partial D_k(x, y)}{\partial y} \cdot f(x, y), \quad k = 1, 2, \ldots \end{aligned} \right.
Espansione della Soluzione con le Funzioni D_k(x, y)
Sostituendo queste funzioni ricorsive nella serie di Taylor, otteniamo l’espansione completa:
\begin{aligned} y(x+h) &= y(x) + h D_1(x, y(x)) + \frac{h^2}{2} D_2(x, y(x)) \\ &\quad + \frac{h^3}{6} D_3(x, y(x)) + \cdots + \frac{h^p}{p!} D_p(x, y(x)) + \mathcal{O}(h^{p+1}). \end{aligned}
Da qui possiamo derivare un metodo numerico che approssima la soluzione della ODE fino all’ordine p:
\left[ \begin{aligned} y_{k+1} &= y_k + h y'_k + \frac{h^2}{2} y''_k + \frac{h^3}{6} y'''_k + \cdots + \frac{h^p}{p!} y^{(p)}_k, \\[1em] y'_k &= D_1(x_k, y_k), \\ y''_k &= D_2(x_k, y_k), \\ y'''_k &= D_3(x_k, y_k), \\ &\vdots \\ y^{(p)}_k &= D_p(x_k, y_k). \end{aligned} \right.
Esempio: Equazione Differenziale con f(x, y) = x + y + xy
Consideriamo la seguente equazione differenziale ordinaria:
\left\{ \begin{aligned} y' &= x + y + xy, \\ y(a) &= y_a. \end{aligned} \right.
Per questa ODE, le funzioni D_k(x, y) si calcolano come segue:
\begin{aligned} D_1(x, y) &= x + y + xy, \\[2em] D_2(x, y) &= \frac{\partial D_1(x, y)}{\partial x} + \frac{\partial D_1(x, y)}{\partial y} \cdot D_1(x, y) \\ &= y + 1 + (x + 1)(xy + x + y), \\[2em] D_3(x, y) &= \frac{\partial D_2(x, y)}{\partial x} + \frac{\partial D_2(x, y)}{\partial y} \cdot D_1(x, y) \\ &= (y + 1)x^3 + (3y + 2)x^2 + (6y + 4)x + 4y + 1, \\[2em] D_4(x, y) &= \frac{\partial D_3(x, y)}{\partial x} + \frac{\partial D_3(x, y)}{\partial y} \cdot D_1(x, y) \\ &= (y + 1)x^4 + (4y + 3)x^3 + (12y + 9)x^2 + (16y + 8)x + 10y + 4, \\[2em] D_5(x, y) &= \frac{\partial D_4(x, y)}{\partial x} + \frac{\partial D_4(x, y)}{\partial y} \cdot D_1(x, y) \\ &= (y + 1)x^5 + (5y + 4)x^4 + (20y + 16)x^3 + (40y + 25)x^2 + (50y + 28)x + 26y + 8, \\[2em] D_6(x, y) &= \frac{\partial D_5(x, y)}{\partial x} + \frac{\partial D_5(x, y)}{\partial y} \cdot D_1(x, y) \\ &= (y + 1)x^6 + (6y + 5)x^5 + (30y + 25)x^4 + (80y + 56)x^3 + (150y + 98)x^2 + (156y + 76)x + 76y + 28. \end{aligned}
Metodo Numerico di Ordine 6
Utilizzando fino alla derivata sesta, possiamo costruire un metodo numerico di ordine 6 per approssimare la soluzione:
\left[ \begin{aligned} y_{k+1} &= y_k + h y'_k + \frac{h^2}{2!} y''_k + \frac{h^3}{3!} y'''_k + \frac{h^4}{4!} y''''_k + \frac{h^5}{5!} y'''''_k + \frac{h^6}{6!} y''''''_k, \\[1em] y'_k &= D_1(x_k, y_k), \\ y''_k &= D_2(x_k, y_k), \\ y'''_k &= D_3(x_k, y_k), \\ y''''_k &= D_4(x_k, y_k), \\ y'''''_k &= D_5(x_k, y_k), \\ y''''''_k &= D_6(x_k, y_k). \end{aligned} \right.
Metodo basato su Taylor con f(x,y) nota migliorato
La costruzione delle funzioni D_k(x,y) richiede la sostituzione di y'(x)=D_1(x,y(x)) che poi viene derivato nel calcolo delle funzioni successive. Se si evita questa sostituzione si ottengono delle funzioni ugualmente calcolabili ma in generale piu semplici.
Consideriamo la sequenza
\begin{aligned} \color{red}y'(x) &= f(x,y(x)) = D_1(x,y(x)) \\ \color{blue}y''(x) &= \dfrac{\mathrm{d}}{\mathrm{d}x}D_1(x,y(x)) \\ &= \dfrac{\partial D_1(x,y(x))}{\partial x} + \dfrac{\partial D_1(x,y(x))}{\partial y} {\color{red}y'(x)} \\ & = D_2(x,y(x),y'(x)) \\ y'''(x) &= \dfrac{\mathrm{d}}{\mathrm{d}x}D_2(x,y(x),y'(x)) \\ &= \dfrac{\partial D_2(x,y(x),y'(x))}{\partial x} + \dfrac{\partial D_2(x,y(x),y'(x))}{\partial y} {\color{red}y'(x)} + \dfrac{\partial D_2(x,y(x),y'(x))}{\partial y'} {\color{blue}y''(x)}\\ & = D_2(x,y(x),y'(x)) \\ &~~\vdots \end{aligned}
questa suggerisce la definizione delle funzioni
D_k(x,y,y',\ldots,y^{(k-1)})
tali che data la soluzione esatta y(x) vale
y^{(k)}(x) = D_k(x,y(x),y'(x),\ldots,y^{(k-1)}(x))
e che soddisfano la ricorrenza:
\left[ \begin{aligned} D_1(x, y) &= f(x, y), \\ D_{k+1}(x, y, y', \ldots, y^{(k)}) &= \frac{\partial D_k(x, y, y', \ldots, y^{(k)})}{\partial x} \\ & + \sum_{j=1}^k \frac{\partial D_k(x, y, y', \ldots, y^{(k-1)})} {\partial y^{(j)}}y^{(j+1)}, \quad k=1,2,\ldots \end{aligned} \right.
Esempio: Equazione Differenziale con f(x, y) = x + y + xy
Consideriamo la seguente equazione differenziale ordinaria:
\left\{ \begin{aligned} y' &= x + y + xy, \\ y(a) &= y_a. \end{aligned} \right.
Per questa ODE, le funzioni D_k(x,\ldots) si calcolano come segue (omettiamo gli argomenti per semplicità)
\begin{aligned} D_1 &= x + y + xy, \\[1em] D_2 &= \frac{\partial D_1}{\partial x} + \frac{\partial D_1}{\partial y} y' \\ &= (x + 1)y' + y, \\[1em] D_3 &= \frac{\partial D_2}{\partial x} + \frac{\partial D_2}{\partial y} y' + \frac{\partial D_2}{\partial y'} y'' \\ &= (x+1)y'' + 2y', \\[1em] D_4 &= \frac{\partial D_3}{\partial x} + \frac{\partial D_3}{\partial y}y' + \frac{\partial D_3}{\partial y'}y''' + \frac{\partial D_3}{\partial y'''}y'''' \\ &= (x+1)y''' + 3y'' \\[1em] D_5 &= \frac{\partial D_4}{\partial x} + \frac{\partial D_4}{\partial y}y' + \frac{\partial D_4}{\partial y'}y'' + \frac{\partial D_4}{\partial y''}y''' + \frac{\partial D_4}{\partial y'''}y'''' \\ &= (x+1)y'''' + 4y''', \\[1em] D_6 &= \frac{\partial D_5}{\partial x} + \frac{\partial D_5}{\partial y}y'+ \frac{\partial D_5}{\partial y'}y''+ \frac{\partial D_5}{\partial y''}y'''+ \frac{\partial D_5}{\partial y'''}y''''+ \frac{\partial D_5}{\partial y''''}y''''' \\ &= (x+1)y'''''+ 5y''''. \end{aligned}
Metodo Numerico di Ordine 6
Ora possiamo costruire un metodo numerico di ordine 6, utilizzando le derivate fino alla sesta:
\left[ \begin{aligned} y_{k+1} &= y_k + h y'_k + \frac{h^2}{2!} y''_k + \frac{h^3}{3!} y'''_k + \frac{h^4}{4!} y''''_k + \frac{h^5}{5!} y'''''_k + \frac{h^6}{6!} y''''''_k, \\[1em] y'_k &= (x_k+1) y_k + x_k, \\ y''_k &= (x_k+1) y'_k + y_k, \\ y'''_k &= (x_k+1) y''_k + 2y'_k, \\ y''''_k &= (x_k+1) y'''_k + 3y''_k, \\ y'''''_k &= (x_k+1) y''''_k + 4y'''_k,, \\ y''''''_k &= (x_k+1) y'''''_k + 5y''''_k,. \end{aligned} \right.
Questo metodo, basato sull’espansione in serie di Taylor, permette di ottenere un’approssimazione di ordine elevato con una complessità computazionale ridotta rispetto alla versione semplice del metodo.