Questo sito utilizza cookies, anche di terze parti, per mostrare pubblicità e servizi in linea con il tuo account. Leggi l'informativa sui cookies.
Username: Password: oppure
C/C++ - integratori numerici... funzionamento errato?? & perfezionamento
Forum - C/C++ - integratori numerici... funzionamento errato?? & perfezionamento

Avatar
pieronaldesi (Normal User)
Newbie


Messaggi: 1
Iscritto: 01/09/2009

Segnala al moderatore
Postato alle 23:19
Martedì, 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();

}

:hail:

Ultima modifica effettuata da pieronaldesi il 01/09/2009 alle 23:31
PM Quote