problems with numerical inegrator!!!
hello
I made that program... bat there are some errors...
how may I correct it??
can you help me???
there are some words written in italian but i think that you can understand anyway. I've removed the file saving of simplettico2 becose of a lack of space!!
thank a lot.
#include <iostream.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;
x+=v*t;
v+=F(x1)*t;
}
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. potenziale sinusoidale \n";
cin>>input;
if (input>3)
{
exit(1);
}
cout <<"\n insert the time interval (I suggest you 0.00001)\n";
cin >>t;
modulo(t);
double o = t*1000000;
int r=0;
cout <<"\n inserirt the total time, please insert a value close to "<<o << "\n";
cin >>T;
numero_volte=T/t;
int d = numero_volte/100;
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();
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();
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();
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();
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();
}
I made that program... bat there are some errors...
how may I correct it??
can you help me???
there are some words written in italian but i think that you can understand anyway. I've removed the file saving of simplettico2 becose of a lack of space!!
thank a lot.
#include <iostream.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;
x+=v*t;
v+=F(x1)*t;
}
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. potenziale sinusoidale \n";
cin>>input;
if (input>3)
{
exit(1);
}
cout <<"\n insert the time interval (I suggest you 0.00001)\n";
cin >>t;
modulo(t);
double o = t*1000000;
int r=0;
cout <<"\n inserirt the total time, please insert a value close to "<<o << "\n";
cin >>T;
numero_volte=T/t;
int d = numero_volte/100;
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();
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();
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();
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();
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();
}
http://www.cplusplus.com/forum/articles/1295
| Don't ask others to debug your broken code without giving a hint what sort of problem they should be searching for. Posting a few hundred lines of code, saying "it doesn't work", will get you ignored. |
ok... you are right!!!
my problem is here...
this is a matrix where i save position speed and energy step by step.
the function eulero transform the old x,v and E in the new after every steps.
then i would write it in a file.
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);
/*
void eulero(double &x, double &v, double t)
{
double x1=x;
x+=v*t;
v+=F(x1)*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();
my problem is here...
this is a matrix where i save position speed and energy step by step.
the function eulero transform the old x,v and E in the new after every steps.
then i would write it in a file.
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);
/*
void eulero(double &x, double &v, double t)
{
double x1=x;
x+=v*t;
v+=F(x1)*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();
Topic archived. No new replies allowed.