http://www.elegio.it/mc2/rk/rk-ode2.html     http://www.elegio.it/mc2/rk/     http://www.elegio.it/calcolatrice/

Istruzioni per l'uso

Mostra le regole per realizzare il file rk-ode2-sistemi.txt

È possibile raccogliere in un unico file diversi problemi distinti in base ad un numero convenzionale scelto a piacere ossia il numero dell'opzione.
Come prima istruzione nel file va definito il vettore opzionivalide[ ] usando una assegnazione, per esempio, di questo tipo:
opzionivalide = [ 28, 5, 21 ];
come scelta questa sarebbe una scelta abbastanza insolita perché di solito si usano numeri interi come 0, 1, 2 etc. Comunque, scelto un dato intero associato ad una opzione bisogna definire una serie di function il cui nome termina con l'intero associato all'opzione. Se indichiamo con # il numero dell'opzione devono venire definite le seguenti function:

• variabili#() == fornisce come risultato il numero delle variabili del sistema. La variabile indipendente deve essere la variabile zeresima e va inclusa nell'elenco.
• bloccodati#() == fornisce come risultato quanti dati, incluse le variabili del sistema e aggiungendo ad esse eventuali variabili ausiliarie, vanno tenute in conto.
• promaschera#() == questa funzione serve per creare la maschera dei dati assegnabili per inizializzare il sistema. Deve servire per inizializzare il vettore delle didascalie ossia titolini[ ] e il vettore dei valori di default ossia valorini[ ].
• sistemazione#() == questa function serve per utilizzare i dati memorizzati nel vettore valorini[ ] ed altri eventuali vettori e variabili ausiliarie.
• inider#(y) == questa function serve per inizializzare il vettore delle variabili di sistema ossia il vettore y[ ] passato in argomento.
• iniprec#(y) == questa funzione serve per far dipendere il criterio di precisione da ogni singola variabile. Se una variabile assume mediamente valori molto più grandi delle altre il suo fattore di precisione deve essere adeguatamente ridotto ( piccolo rispetto ad 1 ) e viceversa ad una variabile con valori piccoli va applicato un fattore di precisione nettamente superiore ad 1.
• sisder#(y,kr,ki) == questa function è quella fondamentale perché specifica come calcolare le derivate prime sulla base dei valori delle variabili del sistema.Il vettore che contiene tutto ossia variabili e loro derivate è il vettore y[ ] mentre l'intero ki specifica dove iniziano le variabili di sistema e l'intero kr specifica dove iniziano le corrispondenti derivate delle variabili del sistema ossia il risultato di questa function.
• confronti#(y) == questa function deve produrre una stringa contenente calcoli fatti a conclusione dell'integrazione. Vengono prelevati valori dal vettore dei risultati finali y[ ] e da altri vettori ausiliari usati per calcolare, se è nota, la soluzione analitica dell'integrazione.
• fapausa#(y) == questa function serve solo in fase di debugging. Normalmente deve restituire una stringa al massimo di lunghezza 1. Fino a quando restituisce una stringa più lunga di 1 il programma fa una stampa di debugging stampando tale stringa a patto che almeno abbia lunghezza maggiore di 3. Dalla prima volta che restituisce una stringa di lunghezza 1 non viene più chiamata.

In una delle function ( ad esempio in inider# ) va definito il vettore instampa[ ] che specifica quali delle variabili vanno stampate e va specificato il valore della variabile tfin che indica l'istante finale della simulazione.
Il tempo deve essere la variabile zeresima del sistema da integrare ( è un vincolo da non dimenticare se non volutamente. L'integrazione termina quando la variabile zeresima assume il valore di tfin e dunque come variabile zeresima potremmo definire una altra variabile ossia una posizione o una velocità... purché questo trucco sia fatto consapevolmente e a ragion veduta )

Un esempio di come realizzare rk-ode2-sistemi.txt

Il file rk-ode2-sistemi.txt potrebbe contenere queste istruzioni

//
//
//  Rungekutta ODE da integrare. 
//
opzionivalide=[0,1,2];
// _____________________________________________
//
//  In cosa consiste l'opzione 0
//
//  Problema del moto piano gravitazionale.
//
function variabili0(){
    return 5;
    }
function bloccodati0(){
    return 10;
    }
//
//  Il sistema di equazioni da integrare...
//
function sisder0(y,kr,ki){
    var c;
    y[kr] = 1;
    y[kr+1] = y[ki+3];
    y[kr+2] = y[ki+4];
    c = - Math.pow( y[ki+1]*y[ki+1] + 
                    y[ki+2]*y[ki+2] + 1.e-50, -1.5 );
    y[kr+3] = c*y[ki+1];
    y[kr+4] = c*y[ki+2];
    }
//
//  Il vettore titolini[ ] contiene le didascalie mentre
//  il vettore valorini[ ] contiene i dati di inizializzazione
//
function promaschera0(){
    titolini[0]="Eccentricità dell'orbita";
    valorini[0]=0.7; 
    titolini[1]="Anomalia eccentrica finale ( all'incirca il tempo finale )";
    valorini[1]="2*PI";
    // alert("Esce da promaschera0()");
    }
//
//  Usa il vettore valorini[ ] per impostare i dati iniziali
//  del sistema di ODE.
//  Per effettuare l'inizializzazione si possono definire variabili
//  e function ausiliarie. In questo caso la variabile ecce e
//  il vettore analitico[ ] costruito tramite la function vero
//
var ecce=0.8;
var analitico=[0];
function vero(ecce,eta){
    // calcola la soluzione analitica per una data eccentricita'
    // ed una data anomalia eccentrica.
    var yv=[0];
    yv[0]=eta-ecce*Math.sin(eta);
    yv[1]=Math.cos(eta)-ecce;
    yv[2]=Math.sin(eta)*Math.sqrt(1-ecce*ecce);
    yv[3]=-Math.sin(eta)/(1-ecce*Math.cos(eta));
    yv[4]=Math.cos(eta)*Math.sqrt(1-ecce*ecce)/(1-ecce*Math.cos(eta));
    return yv;
    }
//
function sistemazione0(){
    if(valorini[0]>=0){if(1>valorini[0]) ecce=valorini[0];};
    tfin=valorini[1];
    analitico=vero(ecce,tfin);
    tfin=analitico[0];
    // alert("Esce da sistemazione0");
    }
//
// Inizializza le derivate e specifica quali vuole stampare
// specificando gli elementi del vettore instampa[ ].
//
function inider0(y){
    y[0]= 0;
    y[1]= 1-ecce; 
    y[2]= 0;
    y[3]= 0;
    y[4]= Math.sqrt((1+ecce)/(1-ecce));
    instampa=[0,1,4]; // ossia stampa le variab. 0, 1 e 4.
    }
//
// Fattore per amplificare o smorzare la precisione richiesta
// per ciascuna variabile.
//
function iniprec0(){
    y[1]=1;
    y[2]=1;
    y[3]=1;
    y[4]=1;
    }
//
//  Calcoli finali di controllo. 
//  Vengono aggiunti dopo i dati dell'ultimo passo.
//
//
function confronti0(y){
    var j,ff="";
    if(1>opzione){
        ff+=xlt+"br/>Soluzione esatta : ";
        for(j=0;instampa.length>j;j++) ff+=" y["+(instampa[j]%n_var)+"]= "+
           analitico[instampa[j]%n_var]+" ;";
        }
    return ff;
    }
//
// Per debuggare eventualmente i primi passi di calcolo.
//
function fapausa0(y){
     // stampa solo stringhe piu' lunghe di tre caratteri.
     // se il risultato è una stringa di 0 o 1 carattere
     // non verra' piu' chiamata. 
     if(totpassi>9) return "0";
     return "y[0]="+y[0]+xlt+"br/>";
     }
// _______________________________________
//
// In cosa consiste l'opzione 1
// 
//
function sisder1(y,kr,ki){     
     y[kr]=1;
     y[kr+1]=-Math.PI*Math.PI*y[ki+2];
     y[kr+2]=y[ki+1];
     }
//
function variabili1(){
    return 3;
    }
function bloccodati1(){
    return 10;
    }
//
function promaschera1(){
    titolini[0]="Tempo di fine simulazione";
    valorini[0]=2;
    // alert("Esce da promaschera1()");
    }
//
function sistemazione1(){
    tfin=valorini[0];
    // alert("Esce da sistemazione1");
    }
//
function inider1(y){
    y[0]= 0;
    y[1]= 0; 
    y[2]= 1;
    instampa=[0,1,2]; // ossia stampa le variab. 0, 1 e 2.
    }
//
function iniprec1(y){
    y[1]=1;
    y[2]=1;
    }
//
function confronti1(y){
    var ff="Con valori di y[0] dispari deve essere y[1]=0 e y[2]=-1 "+
           " se pari allora y[2]=1 ";
    return ff;
    }
//
// Per debuggare eventualmente i primi passi di calcolo.
//
function fapausa1(y){
     // stampa solo stringhe piu' lunghe di tre caratteri.
     // se il risultato è una stringa di 0 o 1 carattere
     // non verra' piu' chiamata. 
     if(totpassi>9) return "0";
     return "y[0]="+y[0]+xlt+"br/>";
     }
//
// ________________________________________
//
//  In cosa consiste l'opzione 2
//
var espofun=3;
function sisder2(y,kr,ki){
     y[kr]=1;
     y[kr+1]=espofun*Math.pow(y[ki],espofun-1);
     }
//
function variabili2(){
    return 2;
    }
function bloccodati2(){
    return 10;
    }
//
function promaschera2(){
    titolini[0]="Tempo di fine simulazione";
    valorini[0]=2;
    titolini[1]="Esponente della soluzione";
    valorini[1]=2;
    // alert("Esce da promaschera2()");
    }
//
function sistemazione2(){
    tfin=valorini[0];
    espofun=valorini[1];
    // alert("Esce da sistemazione2");
    }
//
function inider2(y){
    y[0]= 1;
    y[1]= 1; 
    instampa=[0,1]; // ossia stampa le variab. 0 e 1.
    }
//
function iniprec2(y){
    y[1]=1;
    }
//
function confronti2(y){
    var ff="Dovrebbe ottenere : "+Math.pow(tfin,espofun);
    return ff;
    }
//
// Per debuggare eventualmente i primi passi di calcolo.
//
function fapausa2(y){
     // stampa solo stringhe piu' lunghe di tre caratteri.
     // se il risultato è una stringa di 0 o 1 carattere
     // non verra' piu' chiamata. 
     if(totpassi>9) return "0";
     return "y[0]="+y[0]+xlt+"br/>";
     }
//
// __________________________________________
//

Integratore di sistemi di equazioni differenziali

I sistemi di equazioni differenziali e le function scelte per inizializzarli sono contenuti nel file rk-ode2-sistemi.txt.
Questo integratore di equazioni differenziali alle derivate ordinarie ( ODE == Ordinary Differential Equations ) usa la function setTimeout che evita il blocco del browser e dunque consente di effettuare calcoli di lunga durata.

... ... == Opzione da usare == precisione == minimo numero di passi ammesso == massimo passi consecutivi == massimo passi totali == Quanti riavvii

Risultati

...

...

Annotazioni:0

...

Annotazioni:1

...

Annotazioni:2

...

Annotazioni:3

...

Annotazioni:4

...

Annotazioni:5

...

Annotazioni:6

...

Annotazioni:7

...

Annotazioni:8

...

Annotazioni:9
Il risultato corrisponde al valore del vettore y[ ] al penultimo e all'ultimo passo.

Annotazioni finali

Annotazioni finali

Info sul file dei sistemi da integrare

Inizio del file rk-ode2-sistemi.txt utilizzato

Bibliografia

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 in Fortran, Second Edition ISBN 9780521430647
Qui ne propongo una versione a coefficienti INTERI ( il che elimina qualsiasi possibilità di errore di troncamento nella specifica dell'algoritmo stesso ).

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 y_n rappresenta l'effettivo risultato mentre il vettore y_nt deve essere usato solo a scopo di test; se si discosta eccessivamente da quello di y_n bisogna rieffettuare il calcolo riducendo opportunamente h.