Tredicesima lezione

Sappiamo che moltissimi problemi della fisica possono essere espressi utilizzando equazioni differenziali con derivate ordinarie. Pensiamo all'esempio probabilmente più consueto della dinamica del punto materiale:

che si traduce, in una dimensione, in:

che si può anche scrivere come:

Una generica equazione differenziale alle derivate ordinarie di ordine N si può sempre ridurre, con opportune sostituzioni, ad un sistema di equazioni differenziali del prim'ordine. Nel nostro esempio:

In termini generali, la nostra equazione differenziale di ordine N si può esprimere come:

Dove le funzioni f_i sono note.

Sappiamo inoltre che la soluzione di un'equazione differenziale alle derivate ordinarie non è determinata univocamente dall'equazione stessa, ma richiede opportune condizioni al contorno. Queste (pensiamo al caso dinamico) possono essere espresse come condizioni inizali (x(t₀), v(t₀)), oppure con condizioni iniziali e finali (x(t₀), x(t₁)). Ci occupiamo qui soltanto del caso (più semplice) del problema in cui le condizioni iniziali sono note. Questo permette infatti di ottenere la soluzione facendo evolvere, con metodi opportuni, le condizioni inziali.

Il primo sistema che può venire in mente è quello di sostituire dx con un incremento finito. Questo è il cosiddetto metodo di Eulero. Ad esempio:

oppure, nel caso generale,

Sappiamo però dal caso del calcolo delle derivate che questo metodo presenta alcuni rischi, ed è comunque soggetto ad un'indeterminazione intrinseca dell'ordine del quadrato dell'incremento finito:

Il metodo di Eulero, così com'è, non viene pertanto consigliato per nessun uso pratico.

Esiste un altro metodo, assai popolare nel calcolo scientifico, che permette di approfittare della conoscenza dell'andamento della derivata prima anche in punti intermedi all'incremento considerato. In questo modo la variazione della funzione può essere calcolata anche con precisione di ordine superiore al quadrato dell'incremento.
E' possibile infatti rappresentare esattamente la funzione, utilizzando uno o più valori intermedi della derivata prima, fino ad un certo ordine del suo sviluppo di Taylor. Quest'affermazione non è affatto scontata, anche se molti libri la considerano tale e la descrivono come un sistema appropriato per calcolare una specie di valor medio della derivata sull'incremento finito. Sono indispensabili alcuni conti per accertare la veridicità di tale affermazione: proviamo qui ad abbozzare quelli necessari per ricavare le relazioni che esprimono il metodo di Runge-Kutta al second'ordine, e poi ci fideremo ciecamente dei risultati per il metodo al quart'ordine, che è quello più usato nella pratica.

Supponiamo di cercare la funzione v(t), di cui conosciamo la derivata prima F(t, v(t)). In realtà quello ci interessa è partire dalle condizioni iniziali (v₀, t₀), e ricavare i valori (v₁, v₂, ...) corrispondenti.

Nel problema generale di ordine N vogliamo in particolare riuscire ad esprimere l'incremento di v come una somma di contributi:

Dove ogni contributo si ottiene "raffinando" il contributo precedente in questo modo:

Analizziamo ora il caso particolare del metodo al second'ordine, nel quale vogliamo valutare ed utilizzare solo k₁ e k₂.
Lo sviluppo di Taylor della nostra funzione incognita v(t) vale:

ora, noi conosciamo la derivata prima f della funzione, e dobbiamo fare appello ad un po' di calcolo differenziale per le funzioni di più variabili per esprimere lo sviluppo di Taylor come segue:

Lasciamo ora un attimo da parte questo sviluppo, e concentriamoci sul coefficiente k₂ che, secondo la definizione, e sviluppando in serie la funzione f, può essere espresso in questo modo:

Se ora sostituiamo questi valori di k₁ e k₂ nella nostra espressione iniziale (1), otteniamo:

Confrontando membro a membro questa espressione con lo sviluppo di Taylor di v(t) ottenuto sopra, scopriamo che possiamo scegliere i coefficienti c, a e w in modo da ottenere un incremento che è esattamente uguale allo sviluppo di Taylor al prim'ordine. Questi i coefficienti che realizzano tale uguaglianza:

Il che equivale a porre

Notiamo che possiamo scegliere a piacere il valore di c₂ (che è il punto "intermedio" dell'intervallo considerato). Notiamo inoltre che scegliendo c₂ = 1/2, w₁ si azzera, semplificando in qualche modo i calcoli (ma questa potrebbe anche non essere una scelta ottimale). Il metodo di Runge-Kutta al second'ordine, scegliendo di dividere l'incremento scelto a metà, si riassume quindi in questa semplice espressione:

A questo punto, facciamo finta di aver fatto i conti per il metodo al quart'ordine. Questo è il risultato:

Anche in questo caso, la scelta dei coefficienti presenta una certa arbitrarietà. Qui la scelta è stata diretta ad azzerare il massimo numero di coefficienti possibile: questo riduce il numero di punti intermedi dove occorre conoscere il valore della derivata (che, per le equazioni del second'ordine, deve essere a sua volta calcolata numericamente).

Esiste anche un'altra possibile scelta (coefficienti di Kutta), che offre una migliore resistenza agli errori di troncamento, pur richiedendo un maggior numero di operazioni:

A questo punto, siamo pronti a tradurre in pratica tutto questo?
Quasi: ci resta da capire come scegliere l'incremento iniziale.

Naturalmente si può scegliere un incremento a mano, "piccolo" rispetto all'intervallo di evoluzione richiesto per la funzione. Vediamo in questo come possiamo utilizzare il metodo di Runge-Kutta per scoprire l'evoluzione di un sistema dinamico ad una dimensione (questo programma fa uso di graph_helper.c e graph_helper.h e ). Possiamo modificare il valore della variabile tstep per vedere che effetto fa.

Notiamo inoltre che, utilizzando un incremento costante e non correlato alla variazione della funzione nell'intervallo considerato, sprechiamo energie se la funzione varia poco in quell'intervallo, e corriamo rischi se la funzione varia rapidamente in un piccolo intervallo (esistono funzioni abbastanza patologiche da creare questi problemi). In sostanza nessun algoritmo di integrazione delle equazioni differenziali "serio" è privo di un sistema di modifica dinamica della lunghezza dell'incremento.
Come fare allora per verificare se la variazione della funzione è seguita bene dall'incremento scelto?
Si può provare ad utilizzare un incremento doppio e vedere quanto varia il risultato rispetto a quello ottenuto dalla integrazione con due passi singoli. Infatti possiamo ottenere il valore di

Sia con un passo doppio (il risultato dell'evoluzione è x₁):

che con due passi singoli (il risultato dell'evoluzione è x₂):

La differenza (x₂ - x₁) è una buona stima dell'errore di troncamento che si ottiene col passo scelto, e, come risulta dalle precedenti espressioni, dipende dalla quinta potenza dell'incremento.
Volendo mantenere costante l'errore di troncamento, basta regolare opportunamente il passo.
Vediamo come il nostro esempio di codice cambia utilizzando questa strategia: .

Nell'esempio di programma si può scegliere anche di utilizzare come "forza" quella dell'oscillatore armonico, oppure una forza dipendente dal tempo e dalla velocità come quella espressa dalla Equazione di Duffing:

Questa equazione ha alcune interessanti proprietà, bene evidenziate dalla nostra integrazione numerica. Nonostante l'espressione sia piuttosto complicata, è possibile immaginare un sistema fisico da essa descritto:

Esercizio: Provare ad integrare numericamente una qualche equazione differenziale conosciuta.