[C]Problema integrazione metodo RK4

[torky]

Sto facendo un programma che integri il moto del pendolo doppio,ma ho qualche problema, infatti confrontando con un programma che so per certo funzioni, trovo risultati diversi.Vi volevo chiedere se potete darmi una mano a trovare l'errore visto che a me sembra tutto giusto,l'errore sarà sicuramente concettuale.
Vi posto solamente il ciclo che calcola la traiettoria:

  
for(i=0;i<c->t_max/c->dt;i++){
     t1=p1->theta;
     w1=p1->w;
     t2=p2->theta;
     w2=p2->w;
     tmpt11=w1*c->dt;
     tmpw11=alfa1(t1,t2,w1,w2)*c->dt;
     tmpt12=tmpw11*c->dt;
     tmpw12=alfa1(t1+tmpt11/2,t2,w1+tmpw11/2,w2)*c->dt;
     tmpt13=tmpw12*c->dt;
     tmpw13=alfa1(t1+tmpt12/2,t2,w1+tmpw12/2,w2)*c->dt;
     tmpt14=tmpw13*c->dt;
     tmpw14=alfa1(t1+tmpt13,t2,w1+tmpw13,w2)*c->dt;
     p1->theta+=(tmpt11+2*tmpt12+2*tmpt13+tmpt14)/6.;
     p1->w+=(tmpw11+2*tmpw12+2*tmpw13+tmpw14)/6.;
     tmpt21=w2*c->dt;
     tmpw21=alfa2(t1,t2,w1,w2)*c->dt;
     tmpt22=tmpw21*c->dt;
     tmpw22=alfa2(t1,t2+tmpt21/2,w1,w2+tmpw21/2)*c->dt;
     tmpt23=tmpw22*c->dt;
     tmpw23=alfa2(t1,t2+tmpt22/2,w1,w2+tmpw22/2)*c->dt;
     tmpt24=tmpw23*c->dt;
     tmpw24=alfa2(t1,t2+tmpt23,w1,w2+tmpw23)*c->dt;
     p2->theta+=(tmpt21+2*tmpt22+2*tmpt23+tmpt24)/6.;
     p2->w+=(tmpw21+2*tmpw22+2*tmpw23+tmpw24)/6.;
     t+=c->dt;
     fprintf(stdout,"%lg\t%lg\t%lg\t%lg\t%lg\n",t,p1->theta,p2->theta,p1->w,p2->w);   }

Ciao a tutti!!

[MegaAlchimista]

Allega almeno una spiegazione delle variabili

[torky]

allora:
c->dt è il dt d'integrazione ,c è una struttura che contiene le costanti.
p1 p2 sono due strutture che rappresentano i due pendoli e dentro contengono le variabili:
theta (angolo)
w(velocità angolare);
le funzione alfa,1 e 2, calcolano le accelerazioni angolari
le variabili tmptxy e tmpwxy (x,y € (1,4)) sono i coefficienti del metodo rungekutta quelli che chiamano k1,k2....


Se non capite chiedetemi pure, e grazie per le risposte!

[oregon]

Beh ... ma quali formule hai utilizzato?

[torky]

Originariamente inviato da oregon
Beh ... ma quali formule hai utilizzato?

Quelle che ci sono, non le ho inventate io, per farti un idea guarda qui http://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods

[oregon]

Originariamente inviato da torky
Quelle che ci sono, non le ho inventate io

Questo lo immaginavo ma dato che il problema è tuo, potresti fornire tutti i dettagli del caso, ovvero, come hai provato con il programma che funziona e quali valori ottieni invece tu ... come si sviluppa il tuo programma in relazione alle formule ... ecc ...

[torky]

Originariamente inviato da oregon
Questo lo immaginavo ma dato che il problema è tuo, potresti fornire tutti i dettagli del caso, ovvero, come hai provato con il programma che funziona e quali valori ottieni invece tu ... come si sviluppa il tuo programma in relazione alle formule ... ecc ...

La comparazione rispetto al programma funzionante la faccio confrontando i grafici di theta2 in funzione di theta1,con condizioni iniziali:
theta1=10°=theta2 e
w1=w2=0

quello che il mio programma restituisce è del tutto sballato.

le masse e le lunghezze dei pendoli per ora le ho definite costanti all'interno del programma e quindi non sono importanti.
Per quanto riguarda le formue mi sembrano corrette, faccio quattro piccoli passi intermedi tra una posizione el'altra e poi con queste informazioni calcolo la nuova posizione.Ci deve essere un errore concettuale nell'utilizzo delle formule.Il resto del programma non è necessario postarlo.

Grazie ancora per l'interessamento

[torky]

Se uno un dt due ordini più piccolo le soluzioni vengono identiche però all'altro per avere una soluzione precisa basta un dt dell'ordine 10^-3 a me 10^-5.Questo comporta la crezione di file pesanti e lasciare il computer un minuto a calcolare.
Non capisco perchè questo succeda usiamo tutti e due runge kutta di ordine 4 quindi a parità di dt dovremmo avere la stessa precisione.Per di più ho scritto anche l'integrazione con runge-kutta 2 e verlet e ottengo sempre i medesimi risultati.
Le formule delle accelereazioni le ho controllate e ricontrollate.Non so proprio dove sbattere la capoccia.....