Dimostrazione del Metodo di Runge-Kutta 4

[DavideKern]

Buona giornata. Nel corso di Metodi Computazionali ci è stato illustrato l'algoritmo di Runge-Kutta del quart'ordine per integrare le o.d.e. della forma

\(\displaystyle \begin{cases} \displaystyle \frac{dy}{dx} = f(x, y(x)) \\ \\ y(x_0) = y_0 \end{cases} \)

In particolare, nota una stima \(\displaystyle y_k \) della funzione all'inizio di un intervallo di ampiezza \(\displaystyle h \), il valore successivo \(\displaystyle y_{k+1} \) si ricava iterativamente come

\(\displaystyle y_{k+1} = y_{k} + \frac{1}{6} (k_1 + 2 k_2 + 2 k_3 + k_4) + o(h^5) \)

essendo

\(\displaystyle k_1 = h f(x_k, y_k) \\ k_2 = h f(x_k + h/2, y_k + k_1/2) \\ k_3 = h f(x_k + k/2, y_k + k_2/2) \\ k_4 = h f(x_k + h, y_k + k_3) \)

Il problema è che io non riesco a capire su cosi si basi il metodo. Cercando in rete, il massimo che sono riuscito a capire è che si procede per approssimazioni successive della \(\displaystyle y_{k+1} \) annullando di volta in volta i termini di resto; nel farlo si ricorre a una media pesata delle quattro stime dell'incremento, privilegiando quelle calcolate usando la pendenza "al centro" dell'intervallo (\(\displaystyle k_2, k_3 \), nel nostro caso). Tuttavia mi sembrano considerazioni molto qualitative ed euristiche; neanche sui volumi a mia disposizione si è entrato nel dettaglio tecnico di questo metodo, né nel motivo per cui l'errore locale sia proprio dell'ordine di \(\displaystyle h^5 \). Ho provato a sviluppare in serie di potenze \(\displaystyle y_{k+1} \) e \(\displaystyle y_k \) intorno alla stima \(\displaystyle y_{k + 1/2} \equiv y \left( \frac{x_k + x_{k+1}}{2} \right) \), per poi sottrarre membro a membro (è lo stesso modo con cui si dimostra che il metodo di Eulero-Mid Point/RK2 ha un errore locale di \(\displaystyle o(h^3) \) ). Ho ottenuto

\(\displaystyle y_{k+1} = y_k + h y'_{k + 1/2} + \frac{h^3}{24} y'''_{k + 1/2} + o(h^5) \)

ma penso che questa strategia non porti da nessuna parte.

Qualcun* sa come si può ricavare formalmente questo tipo di approssimazione (ad esempio, ma non solo: con quale criterio si fissano i pesi delle varie stime di incrementi?)? Non sono tanto interessato alla dimostrazione dell'algoritmo di RK in generale (cioè, per un qualunque ordine), però vorrei capire se c'è un modo di dimostrare RK4 in maniera non troppo articolata.

Grazie in anticipo!

[feddy]

Al tempo mi sono trovato bene con questa risposta presa da stack-exchange: https://math.stackexchange.com/question ... tta-method
ci sono diversi link che potrebbero interessarti.

Se posso, su che libri hai guardato? Alcuni testi specializzati su ODE in genere sono molto rigorosi riguardo alla derivazione e all'analisi di metodi RK.

[DavideKern]

I libri a cui mi riferisco sono "Landau, Páez - Computational Problems for Physics" e "Press, Teukolsky, Vetterling, Flannery - Numerical Recipes, The Art of Scientifing Computing", sul quale sembra basarsi il primo (Numerical Recipes è anche presente online. Il metodo di Runge-Kutta è introdotto al capitolo 17.1 qui: http://numerical.recipes/book/book.html). Se conosci libri più specializzati sulla risoluzione numerica di ODE, ogni consiglio è ben accetto!

Grazie mille per il link! Ho capito che \(\displaystyle y_{k+1} \) va sviluppata intorno a \(\displaystyle y_k \) e che vanno uguagliati i coefficienti dello sviluppo a quelli proposti da Runge-Kutta; purtroppo però i dubbi di base restano: per quale motivo è lecito assumere \(\displaystyle y_{k+1} - y_k \) proprio nella forma \(\displaystyle \sum_{i=1}^N \omega_i f(x_k + \nu_i h, y(x_k + \nu_i h)) \)? E perché abbiamo la libertà di scegliere \(\displaystyle \nu_i \equiv 0 \)?

[feddy]

Ora non ho tempo per una risposta ben scritta, ma posso dirti che quei libri sono buoni *solo* se hai gia' chiara la teoria e hai bisogno di vedere velocemente lo schema da usare. Oppure come introduzione al calcolo scientifico attraverso l'implementazione di metodi numerici "semplici". Vado a memoria, in NR mi pare che per RK si limitino per lo piu' a scrivere lo schema numerico per poterlo implementare in C.

In italiano puoi cercare il Quarteroni, Sacco, Saleri - Matematica numerica. Ha un capitolo sulle ODE, quello che cerchi dovrebbe esserci. Altrimenti la bibbia sui metodi per ODE, cioe' "Hairer,Norsett,Wanner - Solving Ordinary Differential Equations I", ma e' in inglese e probabilmente essere ostico da leggere agli inizi.