Runge kutta vs PC
grazie x i suggerimenti, volevo chiarire che
r è un vettore di componenti, considerando il caso di orbite piane, x e y.
r è un vettore di componenti, considerando il caso di orbite piane, x e y.
Risposte
Vuoi un consiglio...
visto che il problema che stai trattando non è particolarmente stiff (sensibile al dato iniziale e ai parametri) e dato che runge kutta è un metodo che richiede ad ogni passo la risoluzione di un sistema non lineare, perchè non usi i classici metodi di eulero? Basterebbe anche eulero in avanti che è esplicito.
Passiamo alla metodologia di implementazione...
Ovviamente, dato che il problema è del 2° ordine, dovrai definire qualcosa del tipo:
z=dr/dt quindi avrai il sistema:
z'=-m/r^2;
r'=z;
A questo punto scrivi le condizioni iniziali come:
z(1)=y'0 (dato iniziale su dr/dt);
r(1)=y0 (dato iniziale su r);
Dato che non mi ricordo a memoria l'espressione di runge kutta del quarto ordine ti faccio un esempio con eulero all'indietro, ma concettualmente devi fare la stessa cosa, semplicemente usando la formula di r-k del 4° ordine.
dove n è il numero di intervalli in cui hai deciso di suddividere il tuo dominio e h è l'ampiezza di ciascun intervallo.
Al posto del comando solve (preso da matlab) puoi usare qualsiasi routine in grado di risolvere un sistema non lineare.
visto che il problema che stai trattando non è particolarmente stiff (sensibile al dato iniziale e ai parametri) e dato che runge kutta è un metodo che richiede ad ogni passo la risoluzione di un sistema non lineare, perchè non usi i classici metodi di eulero? Basterebbe anche eulero in avanti che è esplicito.
Passiamo alla metodologia di implementazione...
Ovviamente, dato che il problema è del 2° ordine, dovrai definire qualcosa del tipo:
z=dr/dt quindi avrai il sistema:
z'=-m/r^2;
r'=z;
A questo punto scrivi le condizioni iniziali come:
z(1)=y'0 (dato iniziale su dr/dt);
r(1)=y0 (dato iniziale su r);
Dato che non mi ricordo a memoria l'espressione di runge kutta del quarto ordine ti faccio un esempio con eulero all'indietro, ma concettualmente devi fare la stessa cosa, semplicemente usando la formula di r-k del 4° ordine.
for i=1:n
solve(' z(i+1)=z(i)+h*(-m/r(i+1)^2)',' r(i+1)=r(i)+h*(z(i+1))','z(i+1),r(i+1)');
enddove n è il numero di intervalli in cui hai deciso di suddividere il tuo dominio e h è l'ampiezza di ciascun intervallo.
Al posto del comando solve (preso da matlab) puoi usare qualsiasi routine in grado di risolvere un sistema non lineare.
non vorrei dire una boiata, ma r è un vettore, quindi bisogna scomporre nelle tre componenti...
Io l'ho presa come l'ha scritta mattoncino. Comunque se il sistema ha geometria sferica il problema non si pone. E' descrivibile in termini della sola r in quanto la dipendenza da fi e theta scompare. Pensa alla legge di gravitazione universale...
Comunque, se dovesse essere usata in modo vettoriale (cosa che non credo proprio), è sufficiente farlo per tutte le componenti.
Comunque, se dovesse essere usata in modo vettoriale (cosa che non credo proprio), è sufficiente farlo per tutte le componenti.
Il fatto è che se la geometria è sferica e l'orbita circolare allora r è cost e varia l'angolo....
e poi chi ci dice che è circolare?
Se consideriamo il campo conservativo gravitazionale ed integraiamo (anche numericamente) le leggi del moto per un corpo di massa m, a seconda del rapporto fra la massa M del corpo che instuara tale campo e la massa m e alle condizioni iniziali del moto si possono avere tre possibili casi: orbita ellittica orbita parabolica orbita iperbolica, e tutto questo solo ed esclusivamente in funzione delle specifiche del problema...
Precisamente si dimostra che se l'energia totale del sistema (cinetica + potenziale), che si conserva sempre perchè se il campo è conservativo l'energia è un integrale primo del moto, ha segno negativo, come per la terra attorno al sole, allora l'orbita è chiusa...etc etc
Pensa alle comete ed hai pianeti....la legge che ne regola il moto "macroscopico" è sempre la stessa, cambiano le condizioni di massa, velocità e posizione inziale...
Una cosa molto bella dei metodi numerici per le eq diff è che in un certo senso le traiettorie che se ne ricavano sono "naturali" nel senso che non sono descrizione geometriche, ma ogni posizione in ogni istante dipemde direttamenre da cosa è successo alla posizione ed all'istante precedente...proprio come in natura
e poi chi ci dice che è circolare?
Se consideriamo il campo conservativo gravitazionale ed integraiamo (anche numericamente) le leggi del moto per un corpo di massa m, a seconda del rapporto fra la massa M del corpo che instuara tale campo e la massa m e alle condizioni iniziali del moto si possono avere tre possibili casi: orbita ellittica orbita parabolica orbita iperbolica, e tutto questo solo ed esclusivamente in funzione delle specifiche del problema...
Precisamente si dimostra che se l'energia totale del sistema (cinetica + potenziale), che si conserva sempre perchè se il campo è conservativo l'energia è un integrale primo del moto, ha segno negativo, come per la terra attorno al sole, allora l'orbita è chiusa...etc etc
Pensa alle comete ed hai pianeti....la legge che ne regola il moto "macroscopico" è sempre la stessa, cambiano le condizioni di massa, velocità e posizione inziale...
Una cosa molto bella dei metodi numerici per le eq diff è che in un certo senso le traiettorie che se ne ricavano sono "naturali" nel senso che non sono descrizione geometriche, ma ogni posizione in ogni istante dipemde direttamenre da cosa è successo alla posizione ed all'istante precedente...proprio come in natura
Sono pienamente daccordo con Giovanni che il moto di un corpo nelo spazio va descritto in funzione di 3 (o 6 se non lo si considera puntiforme) coordinate, tuttavia se il corpo è soggetto alla sola azione di un campo radiale, le componenti della velocità (o meglio, della quantità di moto) in direzioni perpendicolari a quella radiale si conservano.
Comunque, vista l'espressione che ha postato mattoncino, mi sembra che il problema coinvolga la sola variabile r, pertanto credo che una soluzione numerica basata solo sull'integrazione dell'EDO sia quello che stava cercando. Se non fosse così.... Basta che lo specifichi!
E' effettivamente interessante notare che le soluzioni numeriche basate sui metodi numerici esposti dipendano solo da ciò che succede all'istante (o negli istanti, come per i metodi mulistep) precedente. Tuttavia il problema è che il tempo da continuo... diventa discreto! E ad ogni passo di discretizzazione si rischia di perdere in precisione (al limite ottenere una soluzione divergente).
Un quesito per Giovanni: hai giustamente notato che la soluzione numerica di un problema tempo-dipendente si basa solo su ciò che avviene negli istanti precedenti per determinare il comportamento al passo successivo. Consideriamo ora un problema ai limiti (ad esempio un filo casicato con una forza distribuita. E' ancora corretto utilizzare una dipendenza "asimmetrica" del dato successivo da quello precedente? E se invece il problema fosse lo studio della diffusione di inquinante in un fiume con una certa velocità di corrente?
Comunque, vista l'espressione che ha postato mattoncino, mi sembra che il problema coinvolga la sola variabile r, pertanto credo che una soluzione numerica basata solo sull'integrazione dell'EDO sia quello che stava cercando. Se non fosse così.... Basta che lo specifichi!
E' effettivamente interessante notare che le soluzioni numeriche basate sui metodi numerici esposti dipendano solo da ciò che succede all'istante (o negli istanti, come per i metodi mulistep) precedente. Tuttavia il problema è che il tempo da continuo... diventa discreto! E ad ogni passo di discretizzazione si rischia di perdere in precisione (al limite ottenere una soluzione divergente).
Un quesito per Giovanni: hai giustamente notato che la soluzione numerica di un problema tempo-dipendente si basa solo su ciò che avviene negli istanti precedenti per determinare il comportamento al passo successivo. Consideriamo ora un problema ai limiti (ad esempio un filo casicato con una forza distribuita. E' ancora corretto utilizzare una dipendenza "asimmetrica" del dato successivo da quello precedente? E se invece il problema fosse lo studio della diffusione di inquinante in un fiume con una certa velocità di corrente?
Scusate. Posso entrare anch'io nel discorso?
Mi sembra che divaghiate: di ogni problema ne fate
argomento di discussione metodologica e di matematica
astratta.
Il problema posto da mattoncino e' concreto:
determinare il moto di un satellite attorno alla terra.
Allora prima concentriamoci sul problema.
mattoncino butta li' un'equazione di moto generica,
ma per rappresentare il movimento (supposto che avvenga in un piano
che comprende centro della terra e centrom del satellite)
occorre essere piu' precisi.
Intanto occorre sapere dove si trova il satellite al tempo t=0
(supponiamo che disti r0 dal centro della terra, assunto come
origine di coordinate polari, con angolo theta=0). Poi la sua velocita'
iniziale (che supponiamo perpendicolare all'asse 0, e di valore V)
Adesso possiamo scrivere le equazioni che governano il moto, e
precisamente le espressioni di accelerazione trasversale (at) e
radiale (ar): (la faccio breve e do le espressioni)
at = 2*r1*theta1 + r*theta2 = 0 (perche' il moto e' centrale)
ar = r2 - r*thete1^2 = K/r^2 (attrazione newtoniana)
(r1 ed r2 sono rispettivamente la derivata prima e seconda di r
rispetto al tempo, e cosi per theta).
Adesso potete applicare i metodi di soluzione che volete,
per ricavare la traettoria del satellite (in forma polare).
Se non lo fate voi, posso provarci io.
Mi sembra che divaghiate: di ogni problema ne fate
argomento di discussione metodologica e di matematica
astratta.
Il problema posto da mattoncino e' concreto:
determinare il moto di un satellite attorno alla terra.
Allora prima concentriamoci sul problema.
mattoncino butta li' un'equazione di moto generica,
ma per rappresentare il movimento (supposto che avvenga in un piano
che comprende centro della terra e centrom del satellite)
occorre essere piu' precisi.
Intanto occorre sapere dove si trova il satellite al tempo t=0
(supponiamo che disti r0 dal centro della terra, assunto come
origine di coordinate polari, con angolo theta=0). Poi la sua velocita'
iniziale (che supponiamo perpendicolare all'asse 0, e di valore V)
Adesso possiamo scrivere le equazioni che governano il moto, e
precisamente le espressioni di accelerazione trasversale (at) e
radiale (ar): (la faccio breve e do le espressioni)
at = 2*r1*theta1 + r*theta2 = 0 (perche' il moto e' centrale)
ar = r2 - r*thete1^2 = K/r^2 (attrazione newtoniana)
(r1 ed r2 sono rispettivamente la derivata prima e seconda di r
rispetto al tempo, e cosi per theta).
Adesso potete applicare i metodi di soluzione che volete,
per ricavare la traettoria del satellite (in forma polare).
Se non lo fate voi, posso provarci io.
N.D.R.
Nel mio primo intervento (prima di cominciare a divagare...) ho spiegato abbastanza dettagliatamente come risolvere numericamente l'equazione in questione, disinteressandomi del suo preciso significato. Il metodo che ho esposto funziona sia che si parli di corpi in orbita, sia che si parli di prosciutti.
Nel mio primo intervento (prima di cominciare a divagare...) ho spiegato abbastanza dettagliatamente come risolvere numericamente l'equazione in questione, disinteressandomi del suo preciso significato. Il metodo che ho esposto funziona sia che si parli di corpi in orbita, sia che si parli di prosciutti.
il tempo diventa discreto...vero! Però a proposito dei campi conservativi che generano azioni a distanza, come il campo gravitazionale, la velocità con il segnale del campo, se così si può chiamare viaggia alla velocità della luce, quindi il corpo orbitante è "informato" della forza che agisce su di esso e delle sue variazioni in un tempo finito....Detto così è molto in soldoni, ma ho alcuni appunti molto precsisi sul campo gravitazionale in cui questo fenomeno è ben sviscerato...se ti interessa...
Perchè no... sembra interessante!
Intanto che Marco e Giovanni discutono sui massimi sistemi,
fornisco a mattoncino un esempio di soluzione del problema
da lui posto.
E' un esempio di applicazione del metodo Runge-Kutta
realizzato direttamente da MathCad e che utilizza le
equazioni di moto citate la volta scorsa.
La traettoria, dipende dalla velocita' iniziale (parametro
cerchiato in rosso). In questo caso e':
Nauralmente questo e' solo un esempio di cosa si puo' fare
(non e' comunque da usare per lanciare satelliti orbitanti).
Se servono altri dettagli, chiedete.
fornisco a mattoncino un esempio di soluzione del problema
da lui posto.
E' un esempio di applicazione del metodo Runge-Kutta
realizzato direttamente da MathCad e che utilizza le
equazioni di moto citate la volta scorsa.
La traettoria, dipende dalla velocita' iniziale (parametro
cerchiato in rosso). In questo caso e':
Nauralmente questo e' solo un esempio di cosa si puo' fare
(non e' comunque da usare per lanciare satelliti orbitanti).
Se servono altri dettagli, chiedete.
Francamente non ti capisco!
Quale algoritmo usi per risolvere il problema?
(prima hai scritto qualcosa e poi l'hai cancellato)
A me sembra che prima ancora della matematica
occorre avere chiarezza di idee sugli obiettivi
da raggiungere (poi si puo' discutere sui mezzi).
La soluzione da me qui proposta considera prima
una rappresentazione in coordinate polari, che poi
per facilita' di rappresentazione viene convertita
in coordinate cartesiane.
Cosa vuol dire un vettore r funzione di x e y. E poi?
Quale algoritmo usi per risolvere il problema?
(prima hai scritto qualcosa e poi l'hai cancellato)
A me sembra che prima ancora della matematica
occorre avere chiarezza di idee sugli obiettivi
da raggiungere (poi si puo' discutere sui mezzi).
La soluzione da me qui proposta considera prima
una rappresentazione in coordinate polari, che poi
per facilita' di rappresentazione viene convertita
in coordinate cartesiane.
Cosa vuol dire un vettore r funzione di x e y. E poi?
attenzione mattoncino!
ho letto solo superficialmente il tutto, e non entro nel merito, ma mi ha colpito il suggerimento
se per caso tu dovessi provare eulero per il calcolo di un'orbita ellittica di eccentricità non trascurabile, ti consiglio caldamente di fermare il calcolo a metà orbita esatta e di verificare che il punto d'arrivo sia corretto (cosa semplicissima se si parte dal perigeo o dall'apogeo, perchè le coordinate del punto opposto si calcolano in un batter d'occhio).
secondo me Eulero, da perfetto gentiluomo, rischia di doversi dimettere per inefficienza (prevenendo un licenziamento per "giusta causa" vista l'entità dell'errore).
ripeto, mezza orbita, per evitare che errori sistematici si compensino con quelli della seconda mezz'orbita simmetrica.
poi, la seconda verifica alla fine dell'orbita.
sarebbe bello se qualcuno volesse metter in piedi una prova, con diversi passi d'integrazione; io al momento non ce la faccio.
tony
ho letto solo superficialmente il tutto, e non entro nel merito, ma mi ha colpito il suggerimento
quote:
.. perchè non usi i classici metodi di eulero? Basterebbe anche eulero in avanti che è esplicito. [Marco83]
se per caso tu dovessi provare eulero per il calcolo di un'orbita ellittica di eccentricità non trascurabile, ti consiglio caldamente di fermare il calcolo a metà orbita esatta e di verificare che il punto d'arrivo sia corretto (cosa semplicissima se si parte dal perigeo o dall'apogeo, perchè le coordinate del punto opposto si calcolano in un batter d'occhio).
secondo me Eulero, da perfetto gentiluomo, rischia di doversi dimettere per inefficienza (prevenendo un licenziamento per "giusta causa" vista l'entità dell'errore).
ripeto, mezza orbita, per evitare che errori sistematici si compensino con quelli della seconda mezz'orbita simmetrica.
poi, la seconda verifica alla fine dell'orbita.
sarebbe bello se qualcuno volesse metter in piedi una prova, con diversi passi d'integrazione; io al momento non ce la faccio.
tony
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo