To: Studenti_2001/2002 From: "Marco Budinich (tel. 040 6763391)" Subject: Esercitazione 17 gennaio 2002 Cc: Bcc: X-Attachments: :Macintosh G3 HD:472786:Esercizio_n_1.nb: Cari ragazzi, vi allego la file di dimostrazione che abbiamo visto insieme in aula: per far partire il programma basta dare il comando 'Mathematica'. Quando il programma e' partito dovete fare 'Open' della file allegata (Esercizio_n_1.nb) e, quando questa e' aperta, andare nel menu' Kernel>>Evaluation>>Evaluate Initialization. A questo punto e' tutto pronto per lavorare: la descrizione degli esercizi da fare e' dentro la file. Inoltre vi confermo che la lezione del giovedi' si svolgera' dalle 15 alle 16 in aula B. Infine la mia pagina della didattica e': http://www.ts.infn.it/~mbh/Teaching.html . A presto, Marco Budinich To: Studenti_2001/2002 From: "Marco Budinich (tel. 040 6763391)" Subject: Problema con = Cc: Bcc: X-Attachments: Cari ragazzi, giovedi' pomeriggio a lezione ho lasciato in sospeso un problema con l'= in Mathematica. Ecco un esempio chiaro: quadrato = i^2; radice = N[i^(1/2)]; TableForm[Table[{i, quadrato, radice}, {i, 1, 10}]] provate quello che succede e poi ne parliamo a lezione. Cari saluti, Marco Budinich To: Studenti_2001/2002 From: "Marco Budinich (tel. 040 6763391)" Subject: Esercitazione 24 gennaio 2002 Cc: Bcc: X-Attachments: Cari ragazzi, in fondo trovate il programma in C che chiede il raggio e calcola l'area; per eseguirlo salvatelo in una file (per esempio area.c) e poi cosi' lo fate andare: > cc -lm area.c > ./a.out dammi il raggio 1 raggio r = 1.000000 A = 3.141593 > Potete partire da questo esempio per costruire un programma che tabula il seno i cui dati analizzerete con Mathematica. Lo scopo e' fare la differenza fra i dati trovati in C e quelli di Mathematica scritti, per esempio, con 20 cifre significative, per cui, a tutti gli effetti, esatti. Vorrei che generaste e analizzaste i dati in 3 casi: - stampati con 4 cifre dopo la virgola (%.4f) - stampati con 12 cifre dopo la virgola (%.12f) - stampati con 12 cifre dopo la virgola (%.12f) e con x ottenuto dal valore precedente di x sommando 0.1 (per esempio: x = 0.0; for (... { s = sin(x); fprintf.... x = x + 0.1; } ) Come relazione mi mandate per mail (con i nomi dei partecipanti): - il codice c del programma usato, - la copia delle istruzioni di Mathematica usate per l'analisi dei dati, - due righe di commento su cosa avete trovato nel grafico delle differenze nei tre casi elencati. Buon lavoro e a domani, Marco Budinich /* programma che calcola l'area di un cerchio */ #include #include /* istruzioni per il "pre-processing" */ #define pi 3.1415926535897932385 main() { double A, r; /* tutte le variabili dichiarate double */ printf("dammi il raggio "); /* qui chiedo e leggo il raggio NB manca \n */ scanf("%lf", &r); A = pi * r * r; /* il calcolo vero e proprio */ printf("raggio r = %lf A = %lf\n", r, A); /* NB %f stampa i float, %lf i double */ return 0; /* ritornando 0 segnalo a unix che tutto e' ok */ } To: Studenti_2001/2002 From: "Marco Budinich (tel. 040 6763391)" Subject: Esercitazione 31 gennaio 2002 Cc: Bcc: X-Attachments: Cari ragazzi, in fondo trovate due programmi in C che integrano una singola equazione differenziale ordinaria di I ordine con il metodo di Eulero. Il primo lo fa esplicitamente nel for mentre il secondo, un po' piu' elaborato, lo fa con una funzione esterna Eulero(x, y, h). Per il resto i programmi sono identici; in entrambi i casi l'equazione differenziale proposta (y'[x] == x^2 - y[x]/x) e' descritta in f(x, y). Scopo dell'esperienza e' di scrivere nei due programmi il metodo di RungeKutta e successivamente integrare l'equazione. Poi dovete leggere i dati con Mathematica e confrontarli con la soluzione esatta che avrete trovato con DSolve. Come relazione mi mandate per mail (con i nomi dei partecipanti): - il codice c della funzione RungeKutta (meglio se a sua volta essa usa Eulero), - la copia delle istruzioni di Mathematica usate per l'analisi dei dati, - due righe di commento su cosa avete trovato nel grafico delle differenze. - se vi resta tempo potete anche provare a: - studiare come dipende l'errore da h (il passo per l'integrazione numerica) per h molto piccoli, - scrivere una versione di Eulero e RungeKutta che non chiama esplicitamente f(x, y). Infine vedrete che nei programmi in c ci sono alcune cose di cui non abbiamo ancora parlato esplicitamente e precisamente: - una variabile globale (n_chiamate_a_f) e cioe' accessibile da tutte le funzioni; - una variabile carattere (nomefile) che contiene il nome della file inizializzata alla dichiarazione; - una dichiarazione preventiva delle funzioni usate: non e' indispensabile ma e' molto utile (un po' come lavarsi i denti); - un if un po' strano -> if ( (filedati = fopen(nomefile, "w")) == NULL ). Di tutto questo parleremo meglio a lezione ma il funzionamento spero vi sia intuitivo. Buon lavoro e a domani, Marco Budinich /* Versione 1: con passo di Eulero direttamente nel for: la funzione f(x,y) descrive qui una singola equazione differenziale ordinaria di I ordine */ #include /* include le defs. standard dell'input-output */ #include /* include le defs. standard della libreria matematica */ /* Inizio dichiarando esplicitamente tutte le funzioni usate dal programma */ double f(double x, double y); /* ora dichiaro le costanti */ #define x0 0.1 /* il valore iniziale di x */ #define y0 0.5 /* il valore iniziale di y */ #define x1 1.5 /* il valore finale di x */ /* ora dichiaro le variabili globali */ long int n_chiamate_a_f = 0; /* ora il main */ int main(void) { int i; /* l'indice per le iterazioni */ double x, y, h; /* le variabili dell'equazione differenziale */ FILE * filedati; /* la variabile per la file di dati */ char * nomefile = "Eulero_RK.dat"; /* una variabile anche per il nome della file */ /* Iniziamo aprendo il file per la scrittura dei dati */ if ( (filedati = fopen(nomefile, "w")) == NULL ) { printf("non riesco ad aprire la file: %s. Ciao.\n", nomefile); return 1; } /* Inizializzazione delle variabili e stampe */ x = x0; y = y0; h = 0.01; printf ("Inizio, valori iniziali x0 = %lf y0 = %lf h = %lf\n", x, y, h); printf("scrivo nella file: %s\n", nomefile); fprintf (filedati,"%lf %lf\n", x, y); /* Il for per il calcolo */ for (i = 0; x + h <= x1; ++i) { y = y + h * f(x, y); /* Qui si fa un passo di Eulero direttamente nel for */ x = x + h; /* printf ("x = %lf y = %lf\n", x, y); */ fprintf (filedati,"%lf %lf\n", x, y); } /* Fine, ultime cose e via */ printf ("Fine, valori finali x = %lf y = %lf in %d passi chiamate a f %ld\n", x, y, i, n_chiamate_a_f); fclose(filedati); return 0; } /* fine del main, ora le varie funzioni */ double f(double x, double y) /* f definisce l'eq. differenziale */ { /* in questo caso y'[x] == x^2 - y[x]/x */ ++n_chiamate_a_f; return x * x - y / x; } /* Versione 2: con funzione Eulero che a sua volta chiama la funzione f(x,y) che descrive qui una singola equazione differenziale ordinaria di I ordine */ #include /* include le defs. standard dell'input-output */ #include /* include le defs. standard della libreria matematica */ /* Inizio dichiarando esplicitamente tutte le funzioni usate dal programma */ double Eulero(double x, double y, double h); double f(double x, double y); /* ora dichiaro le costanti */ #define x0 0.1 /* il valore iniziale di x */ #define y0 0.5 /* il valore iniziale di y */ #define x1 1.5 /* il valore finale di x */ /* ora dichiaro le variabili globali */ long int n_chiamate_a_f = 0; /* ora il main */ int main(void) { int i; /* l'indice per le iterazioni */ double x, y, h; /* le variabili dell'equazione differenziale */ FILE * filedati; /* la variabile per la file di dati */ char * nomefile = "Eulero_RK.dat"; /* una variabile anche per il nome della file */ /* Iniziamo aprendo il file per la scrittura dei dati */ if ( (filedati = fopen(nomefile, "w")) == NULL ) { printf("non riesco ad aprire la file: %s. Ciao.\n", nomefile); return 1; } /* Inizializzazione delle variabili e stampe */ x = x0; y = y0; h = 0.01; printf ("Inizio, valori iniziali x0 = %lf y0 = %lf h = %lf\n", x, y, h); printf("scrivo nella file: %s\n", nomefile); fprintf (filedati,"%lf %lf\n", x, y); /* Il for per il calcolo */ for (i = 0; x + h <= x1; ++i) { y = y + Eulero(x, y, h); /* Un passo di Eulero */ x = x + h; /* printf ("x = %lf y = %lf\n", x, y); */ fprintf (filedati,"%lf %lf\n", x, y); } /* Fine, ultime cose e via */ printf ("Fine, valori finali x = %lf y = %lf in %d passi chiamate a f %ld\n", x, y, i, n_chiamate_a_f); fclose(filedati); return 0; } /* fine del main, ora le varie funzioni */ double Eulero(double x, double y, double h) /* un 'passo' di Eulero */ { return h * f(x, y); } double f(double x, double y) /* f definisce l'eq. differenziale */ { /* in questo caso y'[x] == x^2 - y[x]/x */ ++n_chiamate_a_f; return x * x - y / x; } To: Studenti_2001/2002 From: "Marco Budinich (tel. 040 6763391)" Subject: Studio errore in funzione di h Cc: Bcc: X-Attachments: Cari ragazzi, oggi a lezione vi ho dato il grafico dello studio di Dy in funzione di h: ecco il programmino in C che calcola i dati necessari. A presto, MB /* Versione 3_h: versione identica alla 3 solo che qui si fa un loop sui varii h */ #include /* include le defs. standard dell'input-output */ #include /* include le defs. standard della libreria matematica */ /* Inizio dichiarando esplicitamente tutte le funzioni usate dal programma */ double RungeKutta(double x, double y, double h, double (*fp)(double x, double y)); double Eulero(double x, double y, double h, double (*fp)(double x, double y)); double f(double x, double y); /* ora dichiaro le costanti (NB tutte maiuscole per distinguerle dalle variabili) */ #define X0 0.1 /* il valore iniziale di x */ #define Y0 0.5 /* il valore iniziale di y */ #define H0 1E-6 /* il valore iniziale di h */ /* ora dichiaro le variabili globali */ long int n_chiamate_a_f = 0; /* ora il main */ int main(void) { int n_passi = 0; /* contatore dei passi fatti */ double x, y, h; /* le variabili dell'equazione differenziale */ FILE * filedati; /* la variabile per la file di dati */ char * nomefile = "Eulero_RK.dat"; /* una variabile anche per il nome della file */ /* Iniziamo aprendo il file per la scrittura dei dati */ filedati = fopen(nomefile, "w"); if ( filedati == NULL ) { printf("non riesco ad aprire la file: %s. Ciao.\n", nomefile); return 1; } /* Inizializzazione delle variabili e stampe */ printf("Inizio, valori iniziali x = %lf y = %lf h = %lf\n", X0, Y0, H0); printf("scrivo nella file: %s\n", nomefile); /* Il for per il calcolo */ for (h = H0; h <= 1E-5; h += 1E-7) { y = Y0 + RungeKutta(X0, Y0, h, f); /* Un passo di Eulero o RungeKutta */ x = X0 + h; fprintf(filedati,"%.18lf %.18lf %.18lf\n", x, y, h); ++n_passi; } /* Fine, ultime cose e via */ printf("Fine, valori finali x = %lf y = %lf in %d passi chiamate a f %ld\n", x, y, n_passi, n_chiamate_a_f); fclose(filedati); return 0; } /* fine del main, ora le varie funzioni */ double Eulero(double x, double y, double h, double (*fp)(double x, double y)) /* un 'passo' di Eulero */ { return h * (*fp)(x, y); } double RungeKutta(double x, double y, double h, double (*fp)(double x, double y)) /* un 'passo' di Runge Kutta */ { return h * (*fp)(x + h / 2, y + Eulero(x, y, h/2, fp)); } double f(double x, double y) /* f definisce l'eq. differenziale */ { /* in questo caso y'[x] == x^2 - y[x]/x */ ++n_chiamate_a_f; return x * x - y / x; } To: Studenti_2001/2002 From: "Marco Budinich (tel. 040 6763391)" Subject: Esercitazione 7 febbraio 2002 Cc: Bcc: X-Attachments: Cari ragazzi, ecco alcuni suggerimenti per l'esercitazione di domani, 7 febbraio 2002. Innanzitutto per risolvere simbolicamente l'equazione con Mathematica vi raccomando di procedere per piccoli passi: ecco una possibile scaletta: - definite l'equazione differenziale, le condizioni iniziali e l'operatore differenziale L - scrivete la soluzione generale dell'omogenea e verificate che soddisfi - Ly == 0 - le condizioni iniziali - scrivete l'integrale particolare di Ly == f con il metodo del Wronskiano - infine scrivete la soluzione generale e verificate soddisfi Ly == f - solo per fare i plot sostituite i valori numerici delle costanti (per es soluzione/.{X0-> 0, V0 -> 1, m -> 1, beta -> 1,....}) Per il programma in C vi allego il main. Ci sono da scrivere le funzioni Eulero e RungeKutta. Alcune osservazioni: - le funzioni ritornano "void" cioe' nulla (in pratica sono l'equivalente delle subroutines Fortran) pero' devono modificare i vettori y e dydx. - cominciate scrivendo Eulero che e' piu' semplice; solo se vi resta tempo fate RungeKutta - per scrivere RungeKutta avrete bisogno di un vettore extra per tenere temporaneamente i valori intermedi di y ad x = x + h / 2 Infine come relazione mi inviate per mail (con i nomi, mi raccomando) - il codice delle sole funzioni Eulero e RungeKutta - la soluzione generale trovata con Mathematica Buon lavoro e a domani, Marco Budinich /* Versione 4: con funzioni Eulero e RungeKutta che a loro volta chiamano la funzione f(x,y,dydx) che descrive un sistema di equazioni differenziali ordinarie di I ordine. La particolare f(x,y,dydx) da usare e' passata dalla chiamata. Il fatto che ora ci siano NEQ equazioni contemporanee (nel nostro caso NEQ = 2) rende piu' complicati sia il main che le funzioni. Per esempio ora la funzione f(x,y,dydx) che descrive le equazioni differenziali non puo' piu' ritornare un singolo valore ma deve modificare il vettore dydx che viene passato alla funzione stessa e la funzione ritorna "void" cioe' nessun valore. Stesso discorso per Eulero e RK. */ #include /* include le defs. standard dell'input-output */ #include /* include le defs. standard della libreria matematica */ /* Inizio dichiarando esplicitamente tutte le funzioni usate dal programma */ void Eulero(double x, int n, double *y, double *dydx, double h, void (*fp)(double x, double *y, double *dydx)); void RungeKutta(double x, int n, double *y, double *dydx, double h, void (*fp)(double x, double *y, double *dydx)); void f(double x, double *y, double *dydx); /* ora dichiaro le costanti (NB tutte maiuscole per distinguerle dalle variabili) */ #define NEQ 2 /* il numero di equazioni differenziali del sistema */ #define T0 0.0 /* il valore iniziale del tempo (variabile indipendente) */ #define X0 0.0 /* il valore iniziale della posizione */ #define V0 2.5 /* il valore iniziale della velocita' */ #define T1 1.5 /* il valore finale del tempo (variabile indipendente) */ #define M 1.0 /* il valore della massa */ #define K 99.0 /* il valore della costante elastica della molla */ #define BETA 0.9 /* il valore del coefficiente di attrito */ /* ora dichiaro le variabili globali */ long int n_chiamate_a_f = 0; /* ora il main */ int main(void) { int n_passi, i; /* contatore dei passi fatti */ double x, y[NEQ], dydx[NEQ]; /* le variabili dell'equazione differenziale */ double h; /* il passo di integrazione */ FILE * filedati; /* la variabile per la file di dati */ char * nomefile = "Eulero_RK.dat"; /* una variabile anche per il nome della file */ /* Iniziamo aprendo il file per la scrittura dei dati */ filedati = fopen(nomefile, "w"); if ( filedati == NULL ) { printf("non riesco ad aprire la file: %s. Ciao.\n", nomefile); return 1; } /* Inizializzazione delle variabili e stampe */ x = T0; y[0] = X0; y[1] = V0; h = 0.01; printf("Inizio, valori iniziali x = %lf h = %lf\n", x, y, h); for (i = 0; i < NEQ; ++i) printf(" y[%d] = %lf", i, y[i]); printf("\nscrivo nella file: %s\n", nomefile); fprintf(filedati,"%lf", x); for (i = 0; i < NEQ; ++i) fprintf(filedati," %lf", y[i]); fprintf(filedati,"\n"); /* Il for per il calcolo */ for (n_passi = 0; x <= T1; ++n_passi) { Eulero(x, NEQ, y, dydx, h, f); /* Un passo di Eulero o RungeKutta: y modificata dalla funzione */ x = x + h; fprintf(filedati,"%lf", x); for (i = 0; i < NEQ; ++i) fprintf(filedati," %lf", y[i]); fprintf(filedati,"\n"); } /* Fine, ultime cose e via */ fclose(filedati); printf("Fine, valori finali x = %lf h = %lf", x, h); for (i = 0; i < NEQ; ++i) printf(" y[%d] = %lf", i, y[i]); printf(" in %d passi chiamate a f %d\n", n_passi, n_chiamate_a_f); return 0; } /* fine del main, ora le varie funzioni */ void Eulero(double x, int n, double *y, double *dydx, double h, void (*fp)(double x, double *y, double *dydx)) /* un 'passo' di Eulero */ { return; } void RungeKutta(double x, int n, double *y, double *dydx, double h, void (*fp)(double x, double *y, double *dydx)) /* un 'passo' di Runge Kutta */ { return; } void f(double x, double *y, double *dydx) /* f definisce il sistema di eq. differenziali */ { ++n_chiamate_a_f; return; } To: Studenti_2001/2002 From: "Marco Budinich (tel. 040 6763391)" Subject: Esercitazione 14 febbraio 2002 Cc: Bcc: X-Attachments: Cari ragazzi, ecco alcuni suggerimenti per l'esercitazione di domani, 14 febbraio 2002. Iniziate dal C per integrare il sistema con - Eulero e RungeKutta (almeno quello di ordine 2) - forza F * sin(OMEGA * t) - forza impulsiva (richiede un main un po' piu' sofisticato) A questo punto leggete e confrontate i dati con le soluzione esatte trovate da Mathematica e: - ricordatevi di fare Chop[] della soluzione trovata con Mathematica per eliminare eventuali piccoli residui immaginari - provate a fare plot parametrici v vs. x del sistema - provate a plottare separatamente le soluzioni dell'omogenea e l'integrale particolare. Infine solito breve mail di relazione e commenti (con i nomi) Buon lavoro e a domani, Marco Budinich To: Studenti_2001/2002 From: "Marco Budinich (tel. 040 6763391)" Subject: Esercitazione 21 febbraio 2002 Cc: Bcc: X-Attachments: Cari ragazzi, alla fine del mail trovate le files per l'esercitazione di domani, 21 febbraio 2002. Inoltre vi includo uno scheletro di main con dichiarazioni ed include alla fine. Se lo desiderate potreste provare odeint con un oscillatore armonico 'nuovo': per esempio vi propongo di usare una forza ad onda quadra che si puo' scrivere in C: float forza(float t) { if (sin(OMEGA * t) > 0) return F; else return -F; } mentre in Mathematica la stessa funzione si puo' scrivere: F (2 UnitStep[Sin[OMEGA t]] - 1) (Se userete questa forza constaterete che Mathematica non e' in grado di risolvere esattamente il problema per cui niente possibilita' di fare confronti con la soluzione esatta domani... Una soluzione impareremo pero' a trovarla la prossima settimana. Qualcuno ha qualche idea ????). Infine mi manderete il solito breve mail di relazione e commenti (con i nomi) che questa volta conterra' solo il main in C che avrete scritto. Buon lavoro e a domani, Marco Budinich =========scheletro del main=========================================================== /* Versione 6: con routine odeint tratta dal libro Numerical Recipes. *** NB *** 1) qui tutti i conti sono fatti in 'float'; 2) i vettori vanno da 1 a N (invece che da 0 a N-1). 3) necessita delle files --> main.c <- la file da modificare (questa) odeint.c <- funzione "interfaccia" rkqs.c <- " "stepper" rkck.c <- " algoritmo nrutil.c <- funzioni di utilita' (vector, matrix etc.) nrutil.h <- include da mettere all'inizio */ #include /* include le defs. standard dell'input-output */ #include /* include le defs. standard della libreria matematica */ #include "nrutil.h" /* include le dichiarazioni delle funzioni di utilita' di Numerical recipes */ /* Inizio dichiarando esplicitamente tutte le funzioni usate dal programma */ void odeint(float *ystart, int nvar, float x1, float x2, float eps, float h1, float hmin, int *nok, int *nbad, void (*fp)(float, float [], float []), void (*rkqs)(float [], float [], int, float *, float, float, float [], float *, float *, void (*)(float, float [], float []))); void rkqs(float y[], float dydx[], int n, float *x, float htry, float eps, float yscal[], float *hdid, float *hnext, void (*fp)(float, float [], float [])); void f(float x, float *y, float *dydx); float forza(float x); /* ora dichiaro le costanti (NB tutte maiuscole per distinguerle dalle variabili) */ #define NEQ 2 /* il numero di equazioni differenziali del sistema */ #define T0 0.0 /* il valore iniziale del tempo (variabile indipendente) */ #define X0 0.0 /* il valore iniziale della posizione */ #define V0 0.0 /* il valore iniziale della velocita' */ #define T1 7.0 /* il valore finale del tempo (variabile indipendente) */ /* ora dichiaro le variabili globali */ float dxsav, *xp, **yp; /* variabili usate da odeint per memorizzare i risultati intermedi */ int kmax, kount; /* contatori usati da odeint per memorizzare il numero di risultati intermedi */ int n_chiamate_a_f; /* il numero delle chiamate alla funzione f */ /* ora il main */ int main(void) { FILE * filedati; /* la variabile per la file di dati */ char * nomefile = "Eulero_RK.dat"; /* una variabile anche per il nome della file */ /* Iniziamo aprendo il file per la scrittura dei dati */ filedati = fopen(nomefile, "w"); if ( filedati == NULL ) { printf("non riesco ad aprire la file: %s. Ciao.\n", nomefile); return 1; } ................. odeint(.....); /* Fine: stampa dei risultati */ ................. /* Chiusura delle files aperte */ fclose(filedati); /* fine... */ return 0; } #include "odeint.c" #include "rkqs.c" #include "rkck.c" #include "nrutil.c" ====================================================================================== To: Studenti_2001/2002 From: "Marco Budinich (tel. 040 6763391)" Subject: Esercitazione 7 marzo 2002 Cc: Bcc: X-Attachments: Cari ragazzi, ecco varie informazioni per l'esercitazione di domani, 7 marzo 2002. Inizio con una lista di siti piu' o meno interessanti dedicati al problema dei tre corpi e simili: http://www.freemars.org/l5/aboutl5.html http://www.ams.org/new-in-math/cover/orbits1.html http://www.irit.fr/COSI/training/complexity-tutorial/henri-poincare.htm http://www.physics.drexel.edu/~steve/triple.html Su http://astro.u-strasbg.fr/~koppen/body/ThreeBody.html potete trovare un applet che simula il problema dei 3 corpi direttamente dal vostro browser. Per chi volesse usare dei numeri realistici ecco qualche dato: #define G 6.67e-11 #define MTERRA 5.98e24 #define MLUNA 7.36e22 #define RTERRALUNA 3.82e8 Ora c'e' la routine RungeKutta4 nel caso non l'abbiate scritta personalmente void RungeKutta4(double x, int n, double *y, double *dydx, double h, void (*fp)(double x, double *y, double *dydx)) { int i; double h2, y_temp[NEQ], dydx_temp_1[NEQ], dydx_temp_2[NEQ], dydx_temp_3[NEQ], dydx_temp_4[NEQ]; h2 = h / 2; (*fp)(x, y, dydx_temp_1); /* (*fp) calcola gli n valori delle derivate al punto (x, y) */ for (i = 0; i < n; ++i) /* qui si memorizza temporaneamente il primo valore di y(x + h/2) */ y_temp[i] = y[i] + h2 * dydx_temp_1[i]; (*fp)(x + h2, y_temp, dydx_temp_2); /* i valori delle derivate al punto (x + h/2, y_temp) */ for (i = 0; i < n; ++i) /* qui si memorizza temporaneamente il secondo valore di y(x + h/2) */ y_temp[i] = y[i] + h2 * dydx_temp_2[i]; (*fp)(x + h2, y_temp, dydx_temp_3); /* i valori delle derivate al punto (x + h/2, y_temp) */ for (i = 0; i < n; ++i) /* qui si memorizza temporaneamente il valore di y(x + h) */ y_temp[i] = y[i] + h * dydx_temp_3[i]; (*fp)(x + h, y_temp, dydx_temp_4); /* i valori delle derivate al punto (x + h, y_temp) */ for (i = 0; i < n; ++i) /* qui si fanno gli n 'passi' di Runge Kutta 4 veri e propri */ { dydx[i] = (dydx_temp_1[i] + 2 * (dydx_temp_2[i] + dydx_temp_3[i]) + dydx_temp_4[i]) / 6; y[i] = y[i] + h * dydx[i]; } return; } In fondo al mail trovate una file di Mathematica (da salvare come testo e aprire successivamente con Mathematica) che contiene celle per leggere i dati e fare diversi plot compresi i "filmini". Nei commenti ci sono un po' di spiegazioni su come fare. Quando avete finito mandatemi per mail solo le routines di C f e quella che calcola la forza totale agente sul corpo i-esimo. Buon lavoro, Marco Budinich (************** Content-type: application/mathematica ************** Mathematica-Compatible Notebook This notebook can be used with any Mathematica-compatible application, such as Mathematica, MathReader or Publicon. The data for the notebook starts with the line containing stars above. To get the notebook into a Mathematica-compatible application, do one of the following: * Save the data starting with the line of stars above into a file with a name ending in .nb, then open the file inside the application; * Copy the data starting with the line of stars above to the clipboard, then use the Paste menu command inside the application. Data for notebooks contains only printable 7-bit ASCII and can be sent directly in email or through ftp in text mode. Newlines can be CR, LF or CRLF (Unix, Macintosh or MS-DOS style). NOTE: If you modify the data for this notebook not in a Mathematica- compatible application, you must delete the line below containing the word CacheID, otherwise Mathematica-compatible applications may try to use invalid cache data. For more information on notebooks and Mathematica-compatible applications, contact Wolfram Research: web: http://www.wolfram.com email: info@wolfram.com phone: +1-217-398-0700 (U.S.) Notebook reader applications are available free of charge from Wolfram Research. *******************************************************************) (*CacheID: 232*) (*NotebookFileLineBreakTest NotebookFileLineBreakTest*) (*NotebookOptionsPosition[ 13274, 332]*) (*NotebookOutlinePosition[ 14128, 359]*) (* CellTagsIndexPosition[ 14084, 355]*) (*WindowFrame->Normal*) Notebook[{ Cell[CellGroupData[{ Cell["Laboratorio n.8 7 III 2002", "Subtitle"], Cell[CellGroupData[{ Cell["Lettura dal file e primi plot dei dati", "Subsubsection"], Cell[BoxData[{ \(\(dati\ = \ Import["\"];\)\), "\[IndentingNewLine]", \(\(Print["\", Length[dati], "\< righe\>"];\)\), "\[IndentingNewLine]", \(\(\({tempo, x1, \ y1, \ z1, \ x2, \ y2, \ z2, \ x3, \ y3, \ z3, \ vx1, \ vy1, \ vz1, \ vx2, \ vy2, \ vz2, \ vx3, \ vy3, \ vz3}\ \ = \ \ Transpose[dati];\)\(\[IndentingNewLine]\) \)\), "\[IndentingNewLine]", \(\(NumeroCorpi\ = \((\ Length[dati[\([1]\)]] - 1)\)/ 6;\)\), "\[IndentingNewLine]", \(\(tmin\ = \ First[tempo];\)\), "\[IndentingNewLine]", \(\(\(tmax\ = \ Last[tempo];\)\(\[IndentingNewLine]\) \)\), "\[IndentingNewLine]", \(\(pltx1\ = \ ListPlot[Transpose[{tempo, x1}], PlotStyle \[Rule] {PointSize[0.011], RGBColor[1, 0, 0]}, \ DisplayFunction \[Rule] Identity];\)\), "\n", \(\(pltx2\ = \ ListPlot[Transpose[{tempo, x2}], PlotStyle \[Rule] {PointSize[0.011], RGBColor[0, 1, 0]}, DisplayFunction \[Rule] Identity];\)\), "\[IndentingNewLine]", \(\(\(pltx3\ = \ ListPlot[Transpose[{tempo, x3}], PlotStyle \[Rule] {PointSize[0.011], RGBColor[0, 0, 1]}, DisplayFunction \[Rule] Identity];\)\(\n\) \)\), "\n", \(\(Show[pltx1, pltx2, \ pltx3, \ DisplayFunction \[Rule] $DisplayFunction];\)\)}], "Input"] }, Closed]], Cell[CellGroupData[{ Cell["Ora le definizioni che occorrono per i \"filmini\"", "Subsubsection"], Cell[BoxData[{ \( (*\ Iniziamo\ con\ le\ definizioni\ \(necessarie : \[IndentingNewLine]\ \ \ \ rng\)\ \ \ \ \ \ \ \ \ \ \ \ \[Rule] \ il\ range\ dei\ plots, \[IndentingNewLine]\ \ \ \ \ \ mm\ \ \ \ \ \ \ \ \ \ \ \ \ \ \[Rule] \ le\ masse\ dei\ pianeti, \[IndentingNewLine]\ \ \ \ \ \ PlanetColors\ \ \ \ \[Rule] \ i\ colori\ dei\ pianeti, \[IndentingNewLine]\ \ \ \ \ \ PlanetRadii\ \ \ \ \ \[Rule] \ i\ raggi\ dei\ pianeti\ \((radice\ cubica\ della\ massa\ ove\ \ possibile)\), \[IndentingNewLine]\ \ \ \ \ \ PosizioniPianeti\ \[Rule] \ la\ matrice\ con\ le\ coppie\ di\ coordinate\ dei\ pianeti, \ \[IndentingNewLine]\ \ \ \ \ \ PlanetGraphics\ \[Rule] \ la\ funzione\ che\ disegna\ i\ pianeti; \ richiede\ come\ \(argomenti : \[IndentingNewLine]\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \(nc : \ \ \ \ il\ numero\ dei\ corpi\ da\ \ disegnare\)\), \[IndentingNewLine]\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ x : \ \ \ \ \ un\ set\ di\ 2\ nc\ coordinate\ \((disegna\ nel\ \(piano\ !\ \))\), \[IndentingNewLine]\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ color : \ i\ colori, \[IndentingNewLine]\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ radii : \ i\ raggi, \[IndentingNewLine]\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ rng : \ \ \ le\ coordinate\ max\ dei\ \(\(disegni\)\(.\)\)\ *) \ \[IndentingNewLine]\[IndentingNewLine]\(mm\ = \ {100, 100, 0.01};\)\), "\[IndentingNewLine]", \(\(rng\ = \ 3;\)\), "\[IndentingNewLine]", \(\(PlanetColors\ = \ Table[Apply[RGBColor[#1, #2, #3] &, RotateRight[{1, 0, 0}, ind]], {ind, 0, NumeroCorpi - 1}];\)\), "\[IndentingNewLine]", \(\(PlanetRadii\ = \ Table[Max[\(\(\@\(mm[\([ind]\)]\)\%3\) rng\)\/100, 0.05], {ind, 1, NumeroCorpi}];\)\), "\[IndentingNewLine]", \(\(\(PosizioniPianeti\ = \ Transpose[{x1, y1, x2, y2, x3, y3}];\)\(\[IndentingNewLine]\) \)\), "\[IndentingNewLine]", \(\(ClearAll[PlanetGraphics];\)\), "\[IndentingNewLine]", \(PlanetGraphics[nc_, \ x_, \ color_, \ radii_, \ rng_]\ := \[IndentingNewLine]Graphics[ Table[{color[\([ind]\)], Disk[{x[\([\((ind - 1)\)*\ 2\ + \ 1]\)], \ x[\([\((ind - 1)\)*\ 2\ + \ 2]\)]}, radii[\([ind]\)]]}, {ind, 1, nc}], \[IndentingNewLine]Frame \[Rule] True, AspectRatio \[Rule] 1, \ PlotRange \[Rule] {{\(-rng\), rng}, {\(-rng\), rng}}]\), "\[IndentingNewLine]", \(\)}], "Input"] }, Closed]], Cell[CellGroupData[{ Cell["Ed infine varie celle per fare plots e \"filmini\"", "Subsubsection"], Cell[BoxData[{ RowBox[{ StyleBox[\( (*\ Per\ iniziare\ questa\ cella\ mostra\ le\ posizioni\ iniziali\ e\ \ finali\ dei\ pianeti\ *) \), FontSize->12], StyleBox["\[IndentingNewLine]", FontSize->14], "\[IndentingNewLine]", \(Show[ PlanetGraphics[NumeroCorpi, First[PosizioniPianeti], \ PlanetColors, \ PlanetRadii, \ rng]];\)}], "\[IndentingNewLine]", \(Show[ PlanetGraphics[NumeroCorpi, Last[PosizioniPianeti], \ PlanetColors, \ PlanetRadii, \ rng]];\)}], "Input"], Cell[BoxData[{ RowBox[{ RowBox[{ StyleBox[\( (*\ Cella\ identica\ alla\ precedente\ solo\ che\ nella\ posizione\ \ finale\ c' e'\ ora\ anche\ la\ traiettoria\ del\ pianeta\ scelto\ *) \), FontSize->12], StyleBox["\[IndentingNewLine]", FontSize->14], "\[IndentingNewLine]", \(Show[ PlanetGraphics[NumeroCorpi, First[PosizioniPianeti], \ PlanetColors, \ PlanetRadii, \ rng]];\)}], "\[IndentingNewLine]"}], "\[IndentingNewLine]", \(PianetaDaTracciare\ \ = \ 3;\), "\[IndentingNewLine]", RowBox[{\(datiTraccia\ = \ Transpose[{\(Transpose[ PosizioniPianeti]\)[\([\((PianetaDaTracciare - 1)\)*2 + 1]\)], \(Transpose[ PosizioniPianeti]\)[\([\((PianetaDaTracciare - 1)\)*2 + 2]\)]}];\), "\[IndentingNewLine]"}], "\[IndentingNewLine]", \(Show[ PlanetGraphics[NumeroCorpi, Last[PosizioniPianeti], \ PlanetColors, \ PlanetRadii, \ rng], Graphics[{PlanetColors[\([PianetaDaTracciare]\)], Line[datiTraccia]}]];\)}], "Input"], Cell[BoxData[ RowBox[{ StyleBox[\( (*\ Qui\ si\ fanno\ tutti\ i\ plot\ dei\ pianeti\ in\ successione . \ Per\ animarli\ dovete\ selezionare\ la\ cella\ con\ tutti\ i\ \ disegni\ e\ dare\ il\ comando\ "\"\ che\ si\ \ trova\ nel\ menu'\ "\"\ *) \), FontSize->12], StyleBox["\[IndentingNewLine]", FontSize->14], "\[IndentingNewLine]", \(For[ind\ = \ 1, \ ind\ \[LessEqual] \ Length[dati], \(++ind\), \[IndentingNewLine]\(Show[\ PlanetGraphics[NumeroCorpi, PosizioniPianeti[\([ind]\)], \ PlanetColors, \ PlanetRadii, \ rng]\ ];\)\[IndentingNewLine]];\)}]], "Input"], Cell[BoxData[{ RowBox[{ StyleBox[\( (*\ Cella\ identica\ alla\ precedente\ solo\ che\ c' e'\ ora\ anche\ la\ traiettoria\ del\ pianeta\ scelto\ *) \), FontSize->12], StyleBox["\[IndentingNewLine]", FontSize->14], "\[IndentingNewLine]", \(PianetaDaTracciare\ = \ 3;\)}], "\[IndentingNewLine]", RowBox[{\(datiTraccia\ = \ Transpose[{\(Transpose[ PosizioniPianeti]\)[\([\((PianetaDaTracciare - 1)\)*2 + 1]\)], \(Transpose[ PosizioniPianeti]\)[\([\((PianetaDaTracciare - 1)\)*2 + 2]\)]}];\), "\[IndentingNewLine]"}], "\[IndentingNewLine]", \(For[ind\ = \ 1, \ ind\ \[LessEqual] \ Length[dati], \(++ind\), \[IndentingNewLine]\(Show[\ PlanetGraphics[NumeroCorpi, PosizioniPianeti[\([ind]\)], \ PlanetColors, \ PlanetRadii, \ rng], Graphics[{PlanetColors[\([PianetaDaTracciare]\)], Line[Take[datiTraccia, ind]]}]\[IndentingNewLine]];\)\[IndentingNewLine]];\)}], \ "Input"] }, Closed]], Cell[CellGroupData[{ Cell[TextData[{ "Qui dentro trovate (", StyleBox["ma solo per i piu' coraggiosi", FontWeight->"Bold"], ") delle celle per applicare una rotazione piana ai dati: in questo modo ci \ si mette nel sistema di riferimento ruotante che ha i primi due pianeti fermi \ (ammesso che ruotino in circolo intorno all'origine)." }], "Text"], Cell[CellGroupData[{ Cell[BoxData[{ \(\(\(Needs["\"]\)\(\[IndentingNewLine]\) \)\), "\[IndentingNewLine]", \(\(datifit\ = \ Transpose[{tempo, x1}];\)\), "\[IndentingNewLine]", \(miofit\ = \ NonlinearFit[datifit, A\ Cos[\[Omega]\ t\ + \[CurlyPhi]\ ], t, {A, \[Omega], \[CurlyPhi]}]\)}], "Input"], Cell[BoxData[ \(2.697361136806192`\ \ Cos[\(\(1.5524830448163751`\)\(\[InvisibleSpace]\)\) + 0.06359583403722949`\ t]\)], "Output"] }, Open ]], Cell[BoxData[{ \(\(Afit\ = \ \(-\(1\/2\)\);\)\), "\[IndentingNewLine]", \(\(\[Omega]fit\ = \ 14.1423;\)\), "\[IndentingNewLine]", \(\(\(\[CurlyPhi]fit\ = \ 0;\)\(\[IndentingNewLine]\) \)\), "\[IndentingNewLine]", \(\(ClearAll[toplot];\)\), "\[IndentingNewLine]", \(\(toplot[t_]\ = \ miofit;\)\ \[IndentingNewLine]\), "\[IndentingNewLine]", \(\(ClearAll[toplot];\)\), "\[IndentingNewLine]", \(\(\(toplot[t_]\ = Afit\ Cos[\[Omega]fit\ t\ + \[CurlyPhi]fit\ \ ];\)\(\[IndentingNewLine]\) \)\), "\[IndentingNewLine]", \(\(pltfit\ = \ Plot[toplot[t], {t, tmin, tmax}, \ DisplayFunction \[Rule] Identity];\)\), "\[IndentingNewLine]", \(\(\(pltdata\ = \ ListPlot[datifit, PlotStyle \[Rule] {PointSize[0.011], RGBColor[1, 0, 0]}, \ DisplayFunction \[Rule] Identity];\)\(\[IndentingNewLine]\) \)\), "\[IndentingNewLine]", \(\(Show[pltdata, pltfit, \ DisplayFunction \[Rule] $DisplayFunction];\)\)}], "Input"], Cell[CellGroupData[{ Cell[BoxData[{\(ClearAll[MatriceRotazione];\), "\[IndentingNewLine]", RowBox[{\(MatriceRotazione[t_]\), " ", "=", " ", RowBox[{\(Sign[\(-Afit\)]\), " ", RowBox[{"(", GridBox[{ {\(Cos[\[Omega]fit\ t\ + \[CurlyPhi]fit]\), \(Sin[\[Omega]fit\ \ t\ + \[CurlyPhi]fit]\)}, {\(-Sin[\[Omega]fit\ t\ + \[CurlyPhi]fit]\), \(Cos[\[Omega]fit\ \ t\ + \[CurlyPhi]fit]\)} }], ")"}]}]}]}], "Input"], Cell[BoxData[ \({{Cos[14.1423`\ t], Sin[14.1423`\ t]}, {\(-Sin[14.1423`\ t]\), Cos[14.1423`\ t]}}\)], "Output"] }, Open ]], Cell[CellGroupData[{ Cell[BoxData[{ \(\(frame\ = \ 9;\)\), "\[IndentingNewLine]", \(MatriceRotazione[frame\ hh] . {x1[\([frame\ + \ 1]\)], y1[\([frame\ + \ 1]\)]}\)}], "Input"], Cell[BoxData[ \({\(-0.4999980788748216`\), 5.9925522948348805`*^-6}\)], "Output"] }, Open ]], Cell[BoxData[{ \(\({x1R, \ y1R}\ = \ Transpose[ Apply[MatriceRotazione[#1] . {#2, #3} &, Transpose[{tempo, x1, y1}], {1}]];\)\), "\[IndentingNewLine]", \(\({x2R, \ y2R}\ = \ Transpose[ Apply[MatriceRotazione[#1] . {#2, #3} &, Transpose[{tempo, x2, y2}], {1}]];\)\), "\[IndentingNewLine]", \(\({x3R, \ y3R}\ = \ Transpose[ Apply[MatriceRotazione[#1] . {#2, #3} &, Transpose[{tempo, x3, y3}], {1}]];\)\)}], "Input"], Cell[BoxData[ \(\(PosizioniPianeti\ = \ Transpose[{x1R, y1R, x2R, y2R, x3R, y3R}];\)\)], "Input"] }, Closed]] }, Open ]] }, FrontEndVersion->"4.1 for Macintosh", ScreenRectangle->{{0, 832}, {0, 604}}, CellGrouping->Manual, WindowSize->{707, 490}, WindowMargins->{{9, Automatic}, {Automatic, 9}}, MacintoshSystemPageSetup->"\<\ 00<0001804P000000`d26_oQon@3:`8g0dL5N`?P0080001804P000000]P2:001 0000I00000400`<30?l00BL?00400@0000000000000006P801T1T00000000000 00000000004000000000000000000000\>" ] (******************************************************************* Cached data follows. If you edit this Notebook file directly, not using Mathematica, you must remove the line containing CacheID at the top of the file. The cache data will then be recreated when you save this file from within Mathematica. *******************************************************************) (*CellTagsOutline CellTagsIndex->{} *) (*CellTagsIndex CellTagsIndex->{} *) (*NotebookFileOutline Notebook[{ Cell[CellGroupData[{ Cell[1727, 52, 48, 0, 61, "Subtitle"], Cell[CellGroupData[{ Cell[1800, 56, 63, 0, 42, "Subsubsection"], Cell[1866, 58, 1406, 27, 249, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[3309, 90, 75, 0, 28, "Subsubsection"], Cell[3387, 92, 2625, 49, 400, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[6049, 146, 75, 0, 28, "Subsubsection"], Cell[6127, 148, 570, 13, 64, "Input"], Cell[6700, 163, 1158, 24, 179, "Input"], Cell[7861, 189, 737, 15, 109, "Input"], Cell[8601, 206, 1144, 24, 193, "Input"] }, Closed]], Cell[CellGroupData[{ Cell[9782, 235, 337, 7, 45, "Text"], Cell[CellGroupData[{ Cell[10144, 246, 337, 6, 67, "Input"], Cell[10484, 254, 146, 3, 25, "Output"] }, Open ]], Cell[10645, 260, 1045, 22, 238, "Input"], Cell[CellGroupData[{ Cell[11715, 286, 448, 8, 52, "Input"], Cell[12166, 296, 124, 2, 25, "Output"] }, Open ]], Cell[CellGroupData[{ Cell[12327, 303, 176, 3, 39, "Input"], Cell[12506, 308, 85, 1, 29, "Output"] }, Open ]], Cell[12606, 312, 526, 12, 53, "Input"], Cell[13135, 326, 111, 2, 25, "Input"] }, Closed]] }, Open ]] } ] *) (******************************************************************* End of Mathematica Notebook file. *******************************************************************) To: Studenti_2001/2002 From: "Marco Budinich (tel. 040 6763391)" Subject: Soluzioni esercitazione 7 marzo 2002 Cc: Bcc: X-Attachments: Cari ragazzi, qui di sotto trovate le funzioni f e accelarazione_totale per la soluzione numerica del problema degli NC corpi (versione Odeint). Inoltre volevo segnalarvi che se usate Odeint e avete tutte le variabili inizialmente a 0 potreste avere qualche problema: una cura che ha funzionato e' modificare la costante TINY in odeint da 1.e-30 a 1.e-10. Infine vi segnalo che parte del materiale dato al corso (mails, files Mathematica, soluzione alcuni esercizi filmini etc.) e' disponibile sul mio sito della didattica http://www.ts.infn.it/~mbh/Teaching.html . Buon lavoro e a presto, Marco Budinich void accelarazione_totale(int i, float * acc_i, float * y); void f(float x, float *y, float *dydx) /* f definisce l'eq. differenziale */ { /* in questo caso il sistema di N corpi */ int i, i3; float acc_i[3], *vy; n_chiamate_a_f++; vy = y + NC * 3; for(i = 1; i <= NC * 3; ++i) dydx[i] = vy[i]; for(i = 1; i <= NC; ++i) { accelarazione_totale(i, acc_i, y); i3 = (i + NC - 1) * 3; dydx[i3 + 1] = acc_i[0]; dydx[i3 + 2] = acc_i[1]; dydx[i3 + 3] = acc_i[2]; } } void accelarazione_totale(int i, float * acc_i, float * y) /* accelarazione_totale da l'accelerazione gravitazionale complessiva sulla massa i-esima */ { int i3, j, j3; float r, r2, fmodulo, deltax, deltay, deltaz; acc_i[0] = acc_i[1] = acc_i[2] = 0.0; i3 = 3 * (i - 1); for (j = 1; j <= NC; ++j) if (j != i) { j3 = 3 * (j - 1); deltax = y[j3 + 1] - y[i3 + 1]; deltay = y[j3 + 2] - y[i3 + 2]; deltaz = y[j3 + 3] - y[i3 + 3]; r2 = deltax * deltax + deltay * deltay + deltaz * deltaz; r = sqrt(r2); fmodulo = m[j - 1] / r2; acc_i[0] += fmodulo * deltax / r; acc_i[1] += fmodulo * deltay / r; acc_i[2] += fmodulo * deltaz / r; } acc_i[0] *= G; acc_i[1] *= G; acc_i[2] *= G; }