Metodi Runge-Kutta
I metodi di Runge-Kutta sono metodi di analisi numerica per l'approssimazione di soluzioni di equazioni differenziali . Sono stati chiamati in onore dei matematici Carl Runge e Martin Wilhelm Kutta , che hanno sviluppato il metodo nel 1901.
Questi metodi si basano sul principio di iterazione , ovvero una prima stima della soluzione viene utilizzata per calcolare una seconda stima più precisa e così via.
Principio generale
Considera il seguente problema:
y′=f(t,y),y(t0)=y0{\ displaystyle y '= f (t, y), \ quad y (t_ {0}) = y_ {0}}
cercheremo di risolvere in un insieme discreto t 0 < t 1 <... < t N . Piuttosto che cercare un metodo diretto, i metodi Runge-Kutta propongono di introdurre i punti intermedi per calcolare per induzione i valori ( t n , y n ) con
{(tnon,io,ynon,io)}1⩽io⩽q{\ displaystyle \ {(t_ {n, i}, y_ {n, i}) \} _ {1 \ leqslant i \ leqslant q}}
tnon,io=tnon+vsio⋅hnon{\ displaystyle t_ {n, i} = t_ {n} + c_ {i} \ cdot h_ {n}}
dove h n = t n + 1 - t n è il passo temporale ec i è nell'intervallo [0; 1] . Per ogni punto intermedio, annotiamo la pendenza corrispondente
pnon,io=f(tnon,io,ynon,io).{\ displaystyle p_ {n, i} = f (t_ {n, i}, y_ {n, i}).}
Quindi, per una soluzione esatta y del problema, abbiamo
y(tnon,io)=y(tnon)+∫tnontnon,iof(t,y(t))dt=y(tnon)+hnon∫0vsiof(tnon+uhnon,y(tnon+uhnon))du,∀io=1,...,q,{\ Displaystyle y (t_ {n, i}) = y (t_ {n}) + \ int _ {t_ {n}} ^ {t_ {n, i}} f (t, y (t)) \ mathrm {d} t = y (t_ {n}) + h_ {n} \ int _ {0} ^ {c_ {i}} f (t_ {n} + uh_ {n}, y (t_ {n} + uh_ {n})) \ mathrm {d} u, \, \ forall i = 1, \ dotsc, q,}
y(tnon+1)=y(tnon)+∫tnontnon+1f(t,y(t))dt=y(tnon)+hnon∫01f(tnon+uhnon,y(tnon+uhnon))du.{\ Displaystyle y (t_ {n + 1}) = y (t_ {n}) + \ int _ {t_ {n}} ^ {t_ {n + 1}} f (t, y (t)) \ mathrm {d} t = y (t_ {n}) + h_ {n} \ int _ {0} ^ {1} f (t_ {n} + uh_ {n}, y (t_ {n} + uh_ {n} )) \ mathrm {d} u.}
Calcoleremo questi integrali con un metodo di quadratura, che possiamo scegliere di essere diverso per due valori distinti di i :
∫0vsiog(u)du≈∑K=1io-1aioKg(vsK),∫01g(u)du≈∑K=1qbKg(vsK),{\ displaystyle \ int _ {0} ^ {c_ {i}} g (u) \ mathrm {d} u \ approx \ sum _ {k = 1} ^ {i-1} a_ {ik} g (c_ { k}), \, \ int _ {0} ^ {1} g (u) \ mathrm {d} u \ approx \ sum _ {k = 1} ^ {q} b_ {k} g (c_ {k} ),}
calcolato qui per g ( u ) = f ( t n + uh n , y ( t n + uh n )) .
Il metodo Runge-Kutta dell'ordine q sarà quindi dato da:
∀io=1,...,q,{tnon,io=tnon+vsiohnon,ynon,io=ynon+hnon∑K=1io-1aioKpnon,Kpnon,io=f(tnon,io,ynon,io){\ displaystyle \ forall i = 1, \ dotsc, q, {\ begin {case} t_ {n, i} & = t_ {n} + c_ {i} h_ {n}, \\ y_ {n, i} & = y_ {n} + h_ {n} \ sum _ {k = 1} ^ {i-1} a_ {ik} p_ {n, k} \\ p_ {n, i} & = f (t_ {n , i}, y_ {n, i}) \ end {case}}}
ynon+1=ynon+hnon∑K=1qbKpnon,K.{\ displaystyle y_ {n + 1} = y_ {n} + h_ {n} \ sum _ {k = 1} ^ {q} b_ {k} p_ {n, k}.}
Il metodo è spesso riassunto dalla tabella dei diversi pesi di quadratura, denominata tabella Butcher :
|
vs1{\ displaystyle c_ {1}}
|
|
vs2{\ displaystyle c_ {2}} |
a21{\ displaystyle a_ {21}}
|
|
vs3{\ displaystyle c_ {3}} |
a31{\ displaystyle a_ {31}} |
a32{\ displaystyle a_ {32}}
|
|
⋮{\ displaystyle \ vdots} |
⋮{\ displaystyle \ vdots} |
|
⋱{\ displaystyle \ ddots}
|
|
vsq{\ displaystyle c_ {q}}
|
aq1{\ displaystyle a_ {q1}}
|
aq2{\ displaystyle a_ {q2}}
|
⋯{\ displaystyle \ cdots}
|
aq,q-1{\ displaystyle a_ {q, q-1}} |
|
|
|
b1{\ displaystyle b_ {1}} |
b2{\ displaystyle b_ {2}} |
⋯{\ displaystyle \ cdots} |
bq-1{\ displaystyle b_ {q-1}} |
bq{\ displaystyle b_ {q}}
|
Il metodo però è coerente .
∀io=2,...,q,∑K=1io-1aioK=vsio{\ displaystyle \ forall i = 2, \ dotsc, q, \ sum _ {k = 1} ^ {i-1} a_ {ik} = c_ {i}}
Esempi
Il metodo Runge-Kutta del primo ordine (RK1)
Questo metodo è equivalente al metodo di Eulero , un semplice metodo per risolvere le equazioni differenziali di 1 ° grado.
Considera il seguente problema:
y′=f(t,y),y(t0)=y0{\ displaystyle y '= f (t, y), \ quad y (t_ {0}) = y_ {0}}Il metodo RK1 utilizza l'equazione
ynon+1=ynon+hf(tnon,ynon){\ displaystyle y_ {n + 1} = y_ {n} + hf \ sinistra (t_ {n}, y_ {n} \ destra)}
dove h è il passaggio dell'iterazione. Il problema è quindi scritto:
|
0{\ displaystyle 0}
|
0{\ displaystyle 0}
|
|
|
1{\ displaystyle 1}
|
Il metodo Runge-Kutta del secondo ordine (RK2)
Il metodo RK2 del punto medio è una composizione del metodo di Eulero :
ynon+1=ynon+hf(tnon+h2,ynon+h2f(tnon,ynon)){\ displaystyle y_ {n + 1} = y_ {n} + hf \ left (t_ {n} + {\ frac {h} {2}}, y_ {n} + {\ frac {h} {2}} f \ sinistra (t_ {n}, y_ {n} \ destra) \ destra)}
dove h è il passaggio dell'iterazione.
Consiste nella stima della derivata a metà della fase di integrazione:
ynon+12=ynon+h2f(tnon,ynon){\ displaystyle y_ {n + {\ frac {1} {2}}} = y_ {n} + {\ frac {h} {2}} f \ sinistra (t_ {n}, y_ {n} \ destra) }
ynon+12′=f(tnon+h2,ynon+12){\ displaystyle y '_ {n + {\ frac {1} {2}}} = f \ left (t_ {n} + {\ frac {h} {2}}, y_ {n + {\ frac {1 } {2}}} \ right)}
e per rifare la fase di completa integrazione a partire da questo preventivo:
ynon+1=ynon+hynon+12′.{\ displaystyle y_ {n + 1} = y_ {n} + hy '_ {n + {\ frac {1} {2}}}.}
Questo schema viene comunemente definito schema correttivo predittivo esplicito .
Questo è il caso particolare di α = 1/2 del metodo più generale:
|
0{\ displaystyle 0} |
0{\ displaystyle 0} |
0{\ displaystyle 0}
|
|
α{\ displaystyle \ alpha}
|
α{\ displaystyle \ alpha}
|
0{\ displaystyle 0}
|
|
|
1-12α{\ displaystyle 1 - {\ tfrac {1} {2 \ alpha}}} |
12α{\ displaystyle {\ tfrac {1} {2 \ alpha}}}
|
Si riconosce così che il metodo di quadratura utilizzato per i tempi intermedi è quello del punto medio .
È un metodo di ordine 2 perché l'errore è dell'ordine di h 3 .
Un altro caso comune è il metodo di Heun, corrispondente al caso α = 1 . Il metodo della quadratura si basa sul metodo del trapezio .
Metodo Runge-Kutta classico del quarto ordine (RK4)
È un caso particolare di uso molto frequente, ha notato RK4.
Considera il seguente problema:
y′=f(t,y),y(t0)=y0{\ displaystyle y '= f (t, y), \ quad y (t_ {0}) = y_ {0}}Il metodo RK4 è dato dall'equazione:
ynon+1=ynon+h6(K1+2K2+2K3+K4){\ displaystyle y_ {n + 1} = y_ {n} + {\ frac {h} {6}} \ sinistra (k_ {1} + 2k_ {2} + 2k_ {3} + k_ {4} \ destra) }
o
K1=f(tnon,ynon){\ displaystyle k_ {1} = f \ sinistra (t_ {n}, y_ {n} \ destra)}
K2=f(tnon+h2,ynon+h2K1){\ displaystyle k_ {2} = f \ sinistra (t_ {n} + {\ frac {h} {2}}, y_ {n} + {\ frac {h} {2}} k_ {1} \ destra) }
K3=f(tnon+h2,ynon+h2K2){\ displaystyle k_ {3} = f \ sinistra (t_ {n} + {\ frac {h} {2}}, y_ {n} + {\ frac {h} {2}} k_ {2} \ destra) }
K4=f(tnon+h,ynon+hK3){\ Displaystyle k_ {4} = f \ sinistra (t_ {n} + h, y_ {n} + hk_ {3} \ right)}
L'idea è che il valore successivo ( y n +1 ) sia approssimato dalla somma del valore corrente ( y n ) e il prodotto della dimensione dell'intervallo ( h ) per la pendenza stimata. La pendenza è ottenuta da una media ponderata delle pendenze:
-
k 1 è la pendenza all'inizio dell'intervallo;
-
k 2 è la pendenza al centro dell'intervallo, utilizzando la pendenza k 1 per calcolare il valore di y nel punto t n + h / 2 mediante il metodo di Eulero ;
-
k 3 è ancora la pendenza nel mezzo dell'intervallo, ma questa volta ottenuta usando la pendenza k 2 per calcolare y ;
-
k 4 è la pendenza alla fine dell'intervallo, con il valore di y calcolato utilizzando k 3 .
Nella media delle quattro piste, maggior peso è dato alle piste nel punto medio.
|
0{\ displaystyle 0} |
0{\ displaystyle 0} |
0{\ displaystyle 0} |
0{\ displaystyle 0} |
0{\ displaystyle 0}
|
|
12{\ displaystyle {\ tfrac {1} {2}}} |
12{\ displaystyle {\ tfrac {1} {2}}} |
0{\ displaystyle 0} |
0{\ displaystyle 0} |
0{\ displaystyle 0}
|
|
12{\ displaystyle {\ tfrac {1} {2}}} |
0{\ displaystyle 0} |
12{\ displaystyle {\ tfrac {1} {2}}} |
0{\ displaystyle 0} |
0{\ displaystyle 0}
|
|
1{\ displaystyle 1}
|
0{\ displaystyle 0}
|
0{\ displaystyle 0}
|
1{\ displaystyle 1}
|
0{\ displaystyle 0}
|
|
|
16{\ displaystyle {\ tfrac {1} {6}}} |
13{\ displaystyle {\ tfrac {1} {3}}} |
13{\ displaystyle {\ tfrac {1} {3}}} |
16{\ displaystyle {\ tfrac {1} {6}}}
|
Il metodo RK4 è un metodo di ordine 4, il che significa che l'errore fatto ad ogni passo è dell'ordine di h 5 , mentre l'errore totale accumulato è dell'ordine di h 4 .
Queste formule sono valide anche per funzioni con valori vettoriali.
Metodo Runge-Kutta del quarto ordine con derivata seconda
Nel caso , possiamo scomporre il problema in un sistema di equazioni:
y″=f(t,y,y′),y(t0)=y0,y′(t0)=y0′{\ displaystyle y '' = f (t, y, y '), \ quad y (t_ {0}) = y_ {0}, \ quad y' (t_ {0}) = y '_ {0}}
{y′=g(t,y,z)z′=S(t,y,z){\ Displaystyle \ left \ {{\ begin {matrix} y '= g (t, y, z) \\ z' = s (t, y, z) \ end {matrix}} \ right.}
con qui y ' = g ( t , y , z ) = z e s = f .
Applicando il metodo RK4 a ciascuna di queste equazioni, quindi semplificando otteniamo:
K1=f(tnon,ynon,ynon′){\ Displaystyle k_ {1} = f \ sinistra (t_ {n}, y_ {n}, y '_ {n} \ destra)}
K2=f(tnon+h2,ynon+h2ynon′,ynon′+h2K1){\ displaystyle k_ {2} = f \ sinistra (t_ {n} + {\ frac {h} {2}}, y_ {n} + {\ frac {h} {2}} y '_ {n}, y '_ {n} + {\ frac {h} {2}} k_ {1} \ right)}
K3=f(tnon+h2,ynon+h2ynon′+h24K1,ynon′+h2K2){\ displaystyle k_ {3} = f \ sinistra (t_ {n} + {\ frac {h} {2}}, y_ {n} + {\ frac {h} {2}} y '_ {n} + {\ frac {h ^ {2}} {4}} k_ {1}, y '_ {n} + {\ frac {h} {2}} k_ {2} \ right)}
K4=f(tnon+h,ynon+hynon′+h22K2,ynon′+hK3){\ displaystyle k_ {4} = f \ sinistra (t_ {n} + h, y_ {n} + hy '_ {n} + {\ frac {h ^ {2}} {2}} k_ {2}, y '_ {n} + HK_ {3} \ right)}
Si deduce y n + 1 e y ' n + 1 grazie a:
ynon+1=ynon+hynon′+h26(K1+K2+K3){\ displaystyle y_ {n + 1} = y_ {n} + hy '_ {n} + {\ frac {h ^ {2}} {6}} \ left (k_ {1} + k_ {2} + k_ {3} \ right)}
ynon+1′=ynon′+h6(K1+2K2+2K3+K4){\ displaystyle y '_ {n + 1} = y' _ {n} + {\ frac {h} {6}} \ left (k_ {1} + 2k_ {2} + 2k_ {3} + k_ {4 } \ giusto)}
Stabilità del metodo
Come metodi a un passo e autoesplicativi, bisogna mettere in discussione la loro stabilità. Possiamo mostrare il seguente risultato:
Supponiamo che f sia k -lipschitziano in y . Sia
I metodi Runge-Kutta sono stabili su [0, T ] , con la costante di stabilità e Λ T , con
α=maxio=1,...,q∑K=1io|aioK|.{\ Displaystyle \ alpha = \ max _ {i = 1, \ dotsc, q} \ sum _ {k = 1} ^ {i} | a_ {ik} |.}
Λ=K∑io=1q|bio|(1+αKhmax+...+(αKhmax)q-1).{\ displaystyle \ Lambda = k \ sum _ {i = 1} ^ {q} | b_ {i} | \ left (1+ \ alpha kh _ {\ max} + \ ldots + (\ alpha kh _ {\ max }) ^ {q-1} \ right).}
Riferimenti
-
[1] Calcoli dettagliati per RK4 con derivata seconda
-
Jean-Pierre Demailly , Analisi numerica ed equazioni differenziali [ dettaglio delle edizioni ]
Articolo correlato
<img src="https://fr.wikipedia.org/wiki/Special:CentralAutoLogin/start?type=1x1" alt="" title="" width="1" height="1" style="border: none; position: absolute;">