Soluzione numerica di sistemi di
equazioni differenziali ordinarie
Sia:
`y = f( y )
Dove f è una funzione del vettore y e `y
rappresenta il vettore delle derivate prime rispetto alla variabile
indipendente solitamente indicata con x o con t.
Si noti che è inutile indicare esplicitamente,
come argomento della funzione f( y ),
la variabile indipendente x perché si può assumere che tale variabile sia la soluzione
dell'equazione differenziale iperbanale:
`x = 1
Dunque basta includere la precedente equazione tra quelle da integrare ed il fatto che
l'algoritmo non incrementi x dell'ampiezza del passo ( ossia che effettuato un passo di ampiezza h,
la variabile x non assuma il valore x+h ) è indice di un sicuro errore di
implementazione dell'algoritmo.
Ovviamente se un algoritmo è dichiarato avere una precisione di ordine O(hp)
esso deve essere capace di risolvere esattamente il seguente sistema:
f(y) = [ 1, p· y(1)p − 1 ]
che, al passo n_esimo, tutti di ampiezza costante h, ha la seguente soluzione:
y(1) = h· n
y(2) = y(1)p
Eulero modificato
Si indica con yn il valore nuovo del vettore
delle variabili y inclusivo della variabile indipendente t.
k1 = h · f( y )
k2 = h · f( y + k1/2 )
yn =
y + k2
|
Questo metodo ha la precisione del secondo ordine ed infatti se poniamo:
y = [ h·n, h2· n2 ]
f(y) = [ 1, 2· y(1) ]
e se proviamo a fare un passo troviamo che:
k1 = h·[ 1, 2·h·n ]
ossia = [ h, 2·h2·n ]
y + k1/2 = [ h· n + h/2 ,
h2· n2 +
h2·n ]
k2 = h · [ 1 , 2·h·n + h ]
ossia = [ h , h2·( 2·n + 1) ]
y + k2 = [ h·n + h ,
h2· n2 +
h2·( 2·n + 1) ]
ossia = [ h·( n + 1) , h2·(
n + 1 )2 ]
come deve essere...
Runge Kutta del quarto ordine
k1 = h · f( y )
k2 = h · f( y + k1/2 )
u = y + k2
k3 = h · f( y + k2/2 )
k4 = h · f( y + k3 )
yn =
y + ( k1 +
2· k2 +
2· k3 +
k4 )/6
|
Notare che u rappresenta già la soluzione essendo il risultato che fornirebbe
il metodo di Eulero modificato e pertanto un criterio per il controllo dell'ampiezza
ottimale del passo potrebbe basarsi sul confronto tra il risultato meno preciso e quello
più preciso yn.
Se la differenza supera un valore ritenuto accettabile si dimezza il passo
e si riprova; se viceversa la concordanza è accettabile, il passo seguente viene
allungato di un certo valore fino a che non si raggiunge l'ampiezza del passo che si
desidera utilizzare stabilmente.
Dato che i coefficienti dell'algoritmo sono numeri razionali è possibile formulare questa
versione a coefficienti interi:
k1 = h · f( y )
k2 = h · f( y + 3 · k1 )
u = y + 6 · k2
k3 = h · f( y + 3 · k2 )
k4 = h · f( y + 6 · k3 )
yn =
y + k1 +
2· k2 +
2· k3 +
k4
|
In questo caso il tempo, ad ogni passo, avanza di 6 · h.
Runge Kutta Nystrom
Citato dall'Abramowitz and Stegun,
Handbook of Mathematical Functions nel cap. 25.5.20 dedicato all'analisi numerica.
Studiato per risolvere soprattutto problemi di dinamica in cui l'accelerazione non
dipende dalla velocità. Esistono però anche versioni in cui l'accelerazione
dipende anche dalla velocità per cui si presume di avere una funzione vettoriale a
dipendente da y e da `y ovvero:
``y = a( y , `y )
Per trattare il tempo, ossia la variabile indipendente, come se fosse una qualsiasi variabile, bisogna
ovviamente porre ``t = 0 ed assegnare come condizioni iniziali t = 0 ; `t = 1.
Si ha:
k1
= h · a( y , `y)
k2
= h · a( y + h · ( `y/2 + k1/8 ) , `y
+ k1/2 )
k3
= h · a( y + h · ( `y/2 + k1/8 ) , `y
+ k2/2 )
k4
= h · a( y + h · ( `y + k3/2 ) , `y
+ k3 )
yn = y + h · ( `y + ( k1
+ k2 + k3 )/6 )
`yn = `y + ( k1 +
2 · k2 + 2 · k3 +
k4 )/6
|
La versione a coefficienti interi di questo algoritmo è la seguente:
k1
= 2 · h · a( y , `y)
k2
= 2 · h · a( y + h · ( 6 · `y + 9 · k1 ) , `y
+ 3 · k1 )
k3
= 2 · h · a( y + h · ( 6 · `y + 9 · k1 ) , `y
+ 3 · k2 )
k4
= 2 · h · a( y + h · ( 12 · `y + 36 · k3 ) , `y
+ 6 · k3 )
yn = y + 12 · h · ( `y + k1
+ k2 + k3 )
`yn = `y + k1 +
2 · k2 + 2 · k3 +
k4
|
La variabile indipendente avanza, in questo caso, di 12 · h.
Il vantaggio di questo algoritmo deriva dal fatto che, quando l'accelerazione
non dipende dalla velocità, il vettore k3
diventa identico al vettore k2 e,
ovviamente, cade la necessità di calcolare `yn;
in sostanza ogni passo richiede tre calcoli dell'accelerazione e non
quattro.
Unico difetto di questo algoritmo è però il fatto che non mi sia nota una
formula per la stima dell'errore.
Risulta successivo (1990) a quello molto noto di Fehlberg reperibile
anche nella Netlib. L'algoritmo di Fehlberg viene citato
come RKF45 e ne ho realizzata
una versione in Fortran 90.
Si tratta di un algoritmo di elevata precisione ( quint'ordine) che fornisce anche
una soluzione di minore precisione (quart'ordine) da usarsi per il controllo del
passo di avanzamento temporale.
Riporto solo la versione a parametri interi:
k1 = h · f( y )
u2 = 16632 ·
k1
k2 = h ·
f( y + 11776 · u2 )
u3 = 6237 ·
( k1 +
3 · k2 )
k3 = h ·
f( y + 11776 · u3 )
u4 = 24948 ·
( k1 −
3 · k2 +
4 · k3 )
k4 = h ·
f( y + 11776 · u4 )
u5 = 1540 ·
( −11 · k1 +
135 · k2 −
140 · k3 +
70 · k4 )
k5 = h ·
f( y + 11776 · u5 )
u6 =
3262 · k1 +
37800 · k2 +
4600 · k3 +
44275 · k4 +
6831 · k5
k6 = h ·
f( y + 8855 · u6 )
u7 =
9361 · k1 +
38500 · k3 +
20125 · k4 +
27648 · k6
yn = y +
10240 · u7
u8 =
39550 · k1 +
148600 · k3 +
94675 · k4 +
7479 · k5 +
96768 · k6
ynt = y +
2530 · u8
|
Ad ogni passo il tempo incrementa di 979292160 · h =
4096·27·5·7·11·23 · h
per cui il valore di h deve essere circa un miliardesimo dell'incremento
che si desidera ottenere.
Il valore di yn rappresenta l'effettivo
risultato mentre il vettore ynt deve
essere usato solo a scopo di test; se si discosta eccessivamente
da quello di yn bisogna rieffettuare
il calcolo riducendo opportunamente h.
Bibliografia:
J.R.Cash, A.H.Karp - ACM Transaction on Mathematical Software,
Vol. 16, No. 3, September 1990, Pages 201-222.
Vedere il
cap. 16.2 di Numerical Recipes
Interpolazione polinomiale cubica
Questo algoritmo è utile quando occorre calcolare il valore del
vettore y ad istanti prefissati, non necessariamente
coincidenti con quelli estremi di un passo di integrazione.
Sia dato il seguente sistema di equazioni differenziali:
`yt = f(yt)
( dove si può assumere che la prima componente di yt sia il tempo t il quale
obbedisce all'equazione differenziale `t = 1 )
e sia stato determinato, con un appropriato metodo
di integrazione:
`yt+h = f(yt+h)
Se 0 ≤ θ ≤ 1 è un parametro adimensionale, allora
l'interpolazione polinomiale cubica (di Hermite) è definita come segue:
yt+θ·h = d1· yt +
d2· `yt +
d3· yt+h +
d4· `yt+h
dove:
d1 = ( θ − 1 )2· ( 2· θ + 1 )
d2 = θ· ( θ − 1 )2 · h
d3 = θ2· ( 3 − 2·θ )
d4 = θ2· ( θ − 1 ) · h
|
|
Notare che solo d4 è sempre non positivo mentre i precedenti tre sono sempre
non negativi.
|