/* leap-frog_HH.c : file sorgente di un programma che consente di osservare l'andamento dell'errore sull'energia per il metodo simplettico di integrazione numerica di tipo leap-frog (o SABA1). Il test viene effettuato controllando la conservazione dell'energia per il modello di Henon e Heiles. */ #include #include /* Definizione dei prototipi delle funzioni */ double energia(double x[2], double y[2], double omega[2]); void kin_flow(double x[2], double y[2], double omega[2], double dt); void pot_flow(double x[2], double y[2], double omega[2], double dt); void saba1(double x[], double y[], double omega[], double c_dt[], double d_dt[]); /* ### */ int main() { int i, numintervalli, i_graph; double x[2], y[2], xiniz[2], yiniz[2], omega[2]; double en_t, en_0, delta_en, errmax_en; double c[2], d[2], Tf, dt, dt_graph, c_dt[2], d_dt[2]; char riga[200]; FILE *ifp, *ofp; printf(" Leggo dal file di input i parametri e le condizioni iniziali\n\n"); /* Inizio delle operazioni di lettura dal file di input */ ifp = fopen("leap-frog_HH.inp" , "r"); /* Leggo "a vuoto" dal file una riga di commenti */ fgets(riga, 200, ifp); /* Leggo dal file una riga di dati */ fgets(riga, 200, ifp); /* Leggo dalla riga le frequenze omega */ sscanf(riga, "%lf %lf", &omega[0], &omega[1]); /* Leggo "a vuoto" dal file una riga di commenti */ fgets(riga, 200, ifp); /* Leggo dal file una riga di dati */ fgets(riga, 200, ifp); /* Leggo dalla riga di dati il valore dell'energia */ sscanf(riga, "%lf", &en_0); /* Leggo "a vuoto" dal file una riga di commenti */ fgets(riga, 200, ifp); /* Leggo dal file una riga di dati */ fgets(riga, 200, ifp); /* Leggo dalla riga di dati i valori iniziali della seconda coppia di x e y */ sscanf(riga, "%lf %lf", &xiniz[1], &yiniz[1]); /* Pongo il valore iniziale della prima coordinata x uguale a 0 */ xiniz[0] = 0; /* Calcolo il primo momento y a partire dall'energia */ delta_en = en_0 - omega[1] * ( pow(xiniz[1],2) + pow(yiniz[1],2) ) / 2 + pow(xiniz[1],3) / 3; yiniz[0] = sqrt( 2 * delta_en / omega[0] ); /* Leggo "a vuoto" dal file una riga di commenti */ fgets(riga, 200, ifp); /* Leggo dal file una riga di dati */ fgets(riga, 200, ifp); /* Leggo dalla riga il tempo finale di integrazione */ sscanf(riga, "%lf", &Tf); /* Leggo "a vuoto" dal file una riga di commenti */ fgets(riga, 200, ifp); /* Leggo dal file una riga di dati */ fgets(riga, 200, ifp); /* Leggo dalla riga il passo di integrazione */ sscanf(riga, "%lf", &dt); /* Leggo "a vuoto" dal file una riga di commenti */ fgets(riga, 200, ifp); /* Leggo dal file una riga di dati */ fgets(riga, 200, ifp); /* Leggo dalla riga l'intervallo di tempo di scansione del grafico */ sscanf(riga, "%lf", &dt_graph); /* Fine delle operazioni di lettura - chiusura del file di input */ fclose(ifp); /* Pongo i valori degli elementi dei vettori c e d uguali alle costanti relative al metodo leap-frog, che e' anche il SABA1, utilizzando la notazione in accordo con l'articolo di Laskar e Robutel del 2001 su Celestial Mechanics and Dynamical Astronomy) */ c[0] = 0.5; d[0] = 1.; /* Pongo i valori delle coordinate canoniche uguali ai corrispondenti valori iniziali letti dal file di input */ for(i=0; i<2; i++) { x[i] = xiniz[i]; y[i] = yiniz[i]; } /* Pongo numintervalli uguale al rapporto tra il tempo finale e il passo di integrazione */ numintervalli = Tf / dt; /* Pongo i_graph uguale al numero di passi che devono intercorrere tra il tracciamento di un punto per il grafico e il successivo */ i_graph = dt_graph / dt; /* Calcolo i valori degli elementi dei vettori c_dt e d_dt in modo che corrispondano ai valori dei passi di integrazione che vengono effettuati (lungo il flusso delle due parti in cui viene divisa la Hamiltoniana) durante l'applicazione del metodo leap-frog */ for(i=0; i<1; i++) { c_dt[i] = c[i] * dt; d_dt[i] = d[i] * dt; } /* Apertura del file di output */ ofp = fopen("leap-frog_HH.out" , "w"); /* Stampo su file esterno il valore dell'istante di tempo iniziale l'errore iniziale sull'energia (entrambi questi valori sono nulli) */ fprintf(ofp, " %le %le\n", 0., 0. ); /* Coerentemente, pongo inizialmente uguale a zero il valore dell'errore massimo sull'energia */ errmax_en = 0.; /* Inizio l'integrazione numerica vera e propria delle equazioni del moto */ printf(" Effettuo l'integrazione numerica delle equazioni del moto\n\n"); for(i=1; i<=numintervalli; i++) { /* Calcolo i valori delle coordinate canoniche al tempo i * dt applicando numintervalli volte il metodo leap-frog */ saba1(x, y, omega, c_dt, d_dt); if(i % i_graph == 0) { /* Calcolo il valore dell'energia corrispondente ai valori assunti dalle coordinate canoniche al tempo i * dt */ en_t = energia(x, y, omega); /* Stampo il valore dell'istante di tempo (i * dt) e la differenza tra l'energia calcolata al tempo i * dt e quella iniziale */ delta_en = en_t - en_0; fprintf(ofp, " %le %le\n", i * dt, delta_en ); /* Aggiorno il valore dell'errore massimo sull'energia */ if(errmax_en < fabs(delta_en)) errmax_en = fabs(delta_en); } } /* Fine delle operazioni di scrittura - chiusura del file di output */ fclose(ofp); /* Stampa finale a proposito della conservazione dell'energia */ printf(" En. iniz. max_t | E(t) - E(0) | max_t | (E(t)-E(0))/E(0) |\n"); printf("%le %le %le\n\n", en_0, errmax_en, fabs(errmax_en/en_0) ); } /* ### */ double energia(double x[2], double y[2], double omega[2]) { int i; double ris = 0; for(i=0; i<2; i++) ris += omega[i] * ( x[i] * x[i] + y[i] * y[i] ) / 2; ris += x[1] * ( x[0] * x[0] - x[1] * x[1] / 3 ); return ris; } /* ### */ void kin_flow(double x[2], double y[2], double omega[2], double dt) { x[0] += omega[0] * y[0] * dt; x[1] += omega[1] * y[1] * dt; } /* ### */ void pot_flow(double x[2], double y[2], double omega[2], double dt) { y[0] -= ( omega[0] * x[0] + 2 * x[0] * x[1] ) * dt; y[1] -= ( omega[1] * x[1] + x[0] * x[0] - x[1] * x[1] ) * dt; } /* ### */ void saba1(double x[2], double y[2], double omega[2], double c_dt[2], double d_dt[2]) { kin_flow(x, y, omega, c_dt[0]); pot_flow(x, y, omega, d_dt[0]); kin_flow(x, y, omega, c_dt[0]); }