pieronaldesi (Normal User)
Newbie
Messaggi: 1
Iscritto: 01/09/2009
|
ciao a tutti
ho provato a fare questo programma si integratori numerici...
mi potete dire se secondo voi funziona bene e perchè non mi stampa l'ultimo file???
come potrei renderli più elegante??
#include <passe_par_tout.h>
double w=1, t, numero_volte, T, m=1, xeu1=1, veu1=1, x2=1, v2=1, xrk2=1, vrk2=1, xrk4=1, vrk4=1, xs1=1, vs1=1, xs2=1, vs2=1;
double Salva_x_v_E [3][101], Salva_x_v_Eeu2 [3][101], Salva_x_v_Erk2 [3][101], Salva_x_v_Erk4 [3][101], Salva_x_v_Es1 [3][101], Salva_x_v_Es2 [3][101];
char *A;
int input;
double F(double x)
{ double Forza;
switch(input)
{
case 1:
Forza = -w*x;
break;
case 2:
Forza = w*x;
break;
case 3:
Forza = -cos(x);
break;int d = numero_volte/100;
}
return Forza;
}
double Fx (double x)
{
double Derivata_Forza;
switch(input)
{
case 1:
Derivata_Forza = -w;
break;
case 2:
Derivata_Forza = w;
break;
case 3:
Derivata_Forza = sin(x);
break;
}
return Derivata_Forza;
}
double energia(double x, double v)
{
double E;
switch(input)
{
case 1:
E=(m*v*v+m*w*w*x*x);
A=(char*)"oscillatore-armonico";
break;
case 2:
E=(m*v*v-m*w*w*x*x);
A=(char*)"oscillatore-iperbolico";
break;
case 3:
E=(m*v*v+m*sin(x));
A=(char*)"potenziale-sinusoidale";
break;
}
return E;
}
void modulo(double &a)
{
a = (a>0)? a :-a;
}
void eulero(double &x, double &v, double t)
{
double x1=x; // x1 diventa l'x del passo precedente!!!!!
x+=v*t; // mi sposto
v+=F(x1)*t; // valuto la velocità
}
void eulero2 (double &x, double &v, double t)
{
double x1=x;
x+=(v + F(x)*t*0.5)*t;
v+=(F(x1) +Fx(x)*v*t*0.5)*t;
}
void runge_kutta2 (double &x, double &v, double t)
{
double K1,K2;
double K1temp, K2temp;
double x1, v1;
x1=x;
v1=v;
K1=v1;
K2=F(x1);
x+=K1*t/2;
v+=K2*t/2;
K1temp=x1+K1*t;
K2temp=v1+K2*t;
K1=K2temp;
K2=F(K1temp);
x+=K1*t/2;
v+=K2*t/2;
}
void runge_kutta4 (double &x, double &v, double t)
{
double K1,K2;
double K1temp, K2temp;
double x1, v1;
x1=x;
v1=v;
K1=v1;
K2=F(x1);
x+=K1*t/6.;
v+=K2*t/6.;
K1temp=x1+K1*t*0.5;
K2temp=v1+K2*t*0.5;
K1=K2temp;
K2=F(K1temp);
x+=K1*t/3.;
v+=K2*t/3.;
K1temp=x1+K1*t*0.5;
K2temp=v1+K2*t*0.5;
K1=K2temp;
K2=F(K1temp);
x+=K1*t/3.;
v+=K2*t/3.;
K1temp=x1+K1*t;
K2temp=v1+K2*t;
K1=K2temp;
K2=F(K1temp);
x+=K1*t/6.;
v+=K2*t/6.;
}
void simplettico1 (double &x, double &v, double t)
{
x+=v*t;
v+=F(x)*t;
}
void simplettico2 (double &x, double &v, double t)
{
x+=v*t/2;
v+=F(x)*t;
x+=v*t/2;
}
int main()
{
cout<<" 1. oscillatore armonico\n 2. oscillatore iperbolico\n 3. pallina puntiforme in un potenziale sinusoidale \n";
cin>>input;
if (input>3)
{
exit(1);
}
cout <<"\n inserire l'intervallo t (valore consigliato circa 0.00001)\n";
cin >>t;
modulo(t);
double o = t*1000000;
int r=0;
cout <<"\n inserire il la durata, ti prego metti un numero dell'ordine di "<<o << "\n";
cin >>T;
numero_volte=T/t;
int d = numero_volte/100;
//eulero primo ordine
Salva_x_v_E [0][0]=xeu1;
Salva_x_v_E [1][0]=veu1;
Salva_x_v_E [2][0]=energia(xeu1,veu1);
for (int i=0; i<numero_volte; i++)
{
eulero(xeu1,veu1,t);
int d = numero_volte/100;
if (i%d==0)
{ r++;
Salva_x_v_E [0][r]=xeu1;
Salva_x_v_E [1][r]=veu1;
Salva_x_v_E [2][r]=energia(xeu1,veu1);
}
}
ofstream Dati ("Dati eulero primo ordine.txt");
Dati <<A<<"\nintegratore eulero primo ordine \n\nl'intervallo di tempo è "<<t<<" \nposizione velocità energia\n";
for (int j=0; j<101; j++)
{Dati<<"\n";
for (int i=0; i<3; i++)
{
Dati<< Salva_x_v_E[j] <<" ";
}
}
Dati.close();
//eulero secondo ordine
Salva_x_v_Eeu2 [0][0]=x2;
Salva_x_v_Eeu2 [1][0]=v2;
Salva_x_v_Eeu2 [2][0]=energia(x2,v2);
for (int i=0; i<numero_volte; i++)
{
eulero2 (x2,v2,t);
int d = numero_volte/100;
if (i%d==0)
{ r++;
Salva_x_v_Eeu2 [0][r]=x2;
Salva_x_v_Eeu2 [1][r]=v2;
Salva_x_v_Eeu2 [2][r]=energia(x2,v2);
}
}
ofstream Dati_eulero_2 ("Dati eulero secondo ordine.txt");
Dati_eulero_2<<A<<"\nintegratore eulero secondo ordine \n\nl'intervallo di tempo è "<<t<<" \nposizione velocità energia\n";
for (int j=0; j<101; j++)
{Dati_eulero_2<<"\n";
for (int i=0; i<3; i++)
{
Dati_eulero_2<< Salva_x_v_Eeu2[j] <<" ";
}
}
Dati_eulero_2.close();
// runge kutta secondo ordine
Salva_x_v_Erk2 [0][0]=xrk2;
Salva_x_v_Erk2 [1][0]=vrk2;
Salva_x_v_Erk2 [2][0]=energia(xrk2,vrk2);
double P, J, N;
for (int i=0; i<numero_volte; i++)
{
runge_kutta2 (xrk2,vrk2,t);
if (i%d==0)
{
r++;
Salva_x_v_Erk2 [0][r]=xrk2;
Salva_x_v_Erk2 [1][r]=vrk2;
Salva_x_v_Erk2 [2][r]=energia(xrk2, vrk2);
}
}
ofstream Dati_rungekutta_2 ("Dati runge kutta secondo ordine.txt");
Dati_rungekutta_2<<A<<"\nintegratore runge kutta secondo ordine\nl'intervallo di tempo è "<<t<<"\nposizione velocità energia\n";
for (int j=0; j<101; j++)
{Dati_rungekutta_2<<"\n";
for (int i=0; i<3; i++)
{
Dati_rungekutta_2<< Salva_x_v_Erk2[j] <<" ";
}
}
Dati_rungekutta_2.close();
// runge kutta quarto ordine
Salva_x_v_Erk4 [0][0]=xrk4;
Salva_x_v_Erk4 [1][0]=vrk4;
Salva_x_v_Erk4 [2][0]=energia(xrk4,vrk4);
for (int i=0; i<numero_volte; i++)
{
runge_kutta4 (xrk4,vrk4,t);
if (i%d==0)
{
r++;
Salva_x_v_Erk4 [0][r]=xrk4;
Salva_x_v_Erk4 [1][r]=vrk4;
Salva_x_v_Erk4 [2][r]=energia(xrk4,vrk4);
}
}
ofstream Dati_rungekutta_4 ("Dati runge kutta quarto ordine.txt");
Dati_rungekutta_4<<A<<"\nintegratore runge kutta quarto ordine \nl'intervallo di tempo è "<<t<<" \nposizione velocità energia\n";
for (int j=0; j<101; j++)
{Dati_rungekutta_4<<"\n";
for (int i=0; i<3; i++)
{
Dati_rungekutta_4<< Salva_x_v_Erk4[j] <<" ";
}
}
Dati_rungekutta_4.close();
// simplettico primo ordine
Salva_x_v_Es1 [0][0]=xs1;
Salva_x_v_Es1 [1][0]=vs1;
Salva_x_v_Es1 [2][0]=energia(xs1,vs1);
for (int i=0; i<numero_volte; i++)
{
simplettico1 (xs1,vs1,t);
if (i%d==0)
{
r++;
Salva_x_v_Es1 [0][r]=xs1;
Salva_x_v_Es1 [1][r]=vs1;
Salva_x_v_Es1 [2][r]=energia(xs1,vs1);
}
}
ofstream Dati_simplettico_1 ("Dati simplettico primo ordine.txt");
Dati_simplettico_1<<A<<"\nintegratore simplettico primo ordine \nl'intervallo di tempo è "<<t<<" \nposizione velocità energia\n";
for (int j=0; j<101; j++)
{Dati_simplettico_1<<"\n";
for (int i=0; i<3; i++)
{
Dati_simplettico_1<< Salva_x_v_Es1[j] <<" ";
}
}
Dati_simplettico_1.close();
// simplettico secondo ordine
Salva_x_v_Es2 [0][0]=xs2;
Salva_x_v_Es2 [1][0]=vs2;
Salva_x_v_Es2 [2][0]=energia(xs2,vs2);
for (int i=0; i<numero_volte; i++)
{
simplettico2 (xs2,vs2,t);
if (i%d==0)
{
r++;
Salva_x_v_Es2 [0][r]=xs2;
Salva_x_v_Es2 [1][r]=vs2;
Salva_x_v_Es2 [2][r]=energia(xs2,vs2);
}
}
ofstream Dati_simplettico_2 ("Dati simplettico secondo ordine.txt");
Dati_simplettico_2<<A<<"\nintegratore simplettico secondo ordine \nl'intervallo di tempo è "<<t<<" \nposizione velocità energia\n";
for (int j=0; j<101; j++)
{Dati_simplettico_2<<"\n";
for (int i=0; i<3; i++)
{
Dati_simplettico_2<< Salva_x_v_Es2[j] <<" ";
}
}
Dati_simplettico_2.close();
}
Ultima modifica effettuata da pieronaldesi il 01/09/2009 alle 23:31 |