[Algoritmi, C] Risoluzione ODE tramite Runge-Kutta 2

[SteezyMenchi]

Salve a tutti. Sono nuovo sul forum di informatica e questo è il mio primo messaggio (si tratta di scientific programming di base in C):
Mi servirebbe aiuto con un integratore di equazioni differenziali: ho messo la scheda negli allegati
Adesso non so se mostrarvi il mio codice qui o direttamente darvi anche quello come file, siccome è abbastanza lungo
Vabbè lo metto in code poi se è necessario uploado il file direttamente.
Ho fatto più o meno quello che mi è stato chiesto (almeno credo): tuttavia, ho implementato il RungeKutta 2nd order non come è scritto nella scheda, ma secondo la notazione del libro Numerical Recipes in C consigliatoci dal professore.
Non so se la mia implementazione sia corretta: ho fatto entrambi i punti 1 e 2(il secondo tramite un file bash e poi plottando con un for loop in gnuplot). è sorto questo problema: se runnate su un compiler online e inserendo i parametri otterrete questo:

N = 0, t0 = 0.010000, x0 = 1.000000, v0 = 0.000000, dt = 0.010000, tmax = 100.000000
t[10000] = 100.0000000000000142 x = 1.5005478528675606 v = -100.0497285075012428

N = 2, t0 = 0.010000, x0 = 1.000000, v0 = 0.000000, dt = 0.010000, tmax = 100.000000
t[10000] = 100.0000000000000142 x = 0.6636376775159283 v = 0.1073845956605055

N = 1, t0 = 0.010000, x0 = 1.000000, v0 = 0.000000, dt = 0.010000, tmax = 100.000000
t[10000] = 100.0000000000000142 x = 0.6030424837371133 v = 0.0696976399374857

N = 3, t0 = 0.010000, x0 = 1.000000, v0 = 0.000000, dt = 0.010000, tmax = 100.000000
t[10000] = 100.0000000000000142 x = 0.7045182928011882 v = 0.1050771817756628

N = 4, t0 = 0.010000, x0 = 1.000000, v0 = 0.000000, dt = 0.010000, tmax = 100.000000
t[10000] = 100.0000000000000142 x = 0.7345036132091616 v = 0.1031850670880924

N = 5, t0 = 0.010000, x0 = 1.000000, v0 = 0.000000, dt = 0.010000, tmax = 100.000000
t[10000] = 100.0000000000000142 x = 0.7577350957939557 v = 0.1015888677781990

N = 6, t0 = 0.010000, x0 = 1.000000, v0 = 0.000000, dt = 0.010000, tmax = 100.000000
t[10000] = 100.0000000000000142 x = 0.7764213062976668 v = 0.1002059937837422

N = 7, t0 = 0.010000, x0 = 1.000000, v0 = 0.000000, dt = 0.010000, tmax = 100.000000
t[10000] = 100.0000000000000142 x = 0.7918693982489068 v = 0.0989846199936164

Fin qui mi sembrava tutto ok. Poi però ho plottato la coordinata $P$ in funzione di $r$ ( io nel codice ho definito $P, D, r$ come $x, v, t$ per abitudine) e ho notato che per tutti gli $n$ da $0 \to 7$ nessuna delle curve interseca l'asse x (ovvero P(che è il mio x) non si annulla mai nel tempo di integrazione).
Tuttavia nel punto dopo chiede di trovare i valori di $R$(ovvero $t$) per cui $P(R)$ si annulla
Se qualcuno potesse gentilmente dirmi dove e cosa sto sbagliando (credo che non ci siano errori nella funzione dove implemento l'algoritmo, forse sbaglio altrove). Ho fatto un po' di debugging con dei printf e mi è sembrato che tutto fosse giusto. Grazie mille in anticipo.
Vorrei specificare che il mio livello di conoscenza del C è davvero basso e si limita alle cose presenti nel programma da me scritto, quindi vi chiedo, nel caso di spiegazioni, di non essere troppo crudeli, e soprattutto di suggerire correzioni all'interno del file.c che ho messo.
Ringrazio in anticipo chi avrà voglia di darmi una mano f1
(P.S.ho iniziato da pochissimo gli algoritmi di integrazione di ODEs)

Codice:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>




/***********STRUCTS**************/


struct Point_t {
  double x;
  double v;
};
typedef struct Point_t Point;

struct Params_t {

  int n;

};
typedef struct Params_t Params;




/*********STRUCT FUNCTIONS DECLARATION***********/

void pSum(Point *res, Point *sum);
void pMul(double factor, Point *res);
void HPMul(double factor, Point *res);
void pFma(double factor, Point *res, Point *p1, Point *p2);
void pCpy(Point *pIn, Point *pOut);
void pSet(double val1, double val2 , Point *pIn);
void printP(Point *pIn);

/*********PROBLEM FUNCTIONS DECLARATION**********/

void drift(Point *pIn, Point *pOut, double t, Params *parameters);
void RK2(Point *pIn, Point *pOut, double t, double dt, Params *parameters);


/***********MAIN**************/

int main(int argc, char *argv[]) {
   
    FILE *fp;
    int i = 0;
    int N = atoi(argv[1]);
    double t0 = atof (argv[2]);
    double x0 = atof(argv[3]);
    double v0 = atof(argv[4]);
    double dt = atof(argv[5]);
    double tmax = atof(argv[6]);
   
    printf("N = %d, t0 = %lf, x0 = %lf, v0 = %lf, dt = %lf, tmax = %lf\n", N, t0, x0, v0, dt, tmax);
   
    double numberofSteps = (int) (tmax / dt -0.5) ;   // approximate to the nearest bigger integer
   
    Point old_xv, new_xv;
    //Assign the initial values
    old_xv.x = x0;
    old_xv.v = v0;
   
   
    Point *pIn, *pOut;
   
    pIn = &old_xv;
    pOut = &new_xv;
   
    Params parameters;
    parameters.n = N;
   
    /*test of drift functions-----------------drift(pIn, pOut, t0, &parameters); works*/
   
    /*Point *null;
    Point zero;
    zero.x = 0;
    zero.v = 0;
    null = &zero;*/
   
    /* Test the first iteration with known outputs RK2(pIn,pOut, t0, dt, &parameters); works*/
   
    for (i = 0; i <= numberofSteps; ++i) {
       
        Point *tmp;
        double t = t0 + dt*i;
        //printf("t[%d] = %.16lf x[%d] = %.16lf v[%d] = %.16lf\n", i, t, i, pIn->x, i, pIn->v);
        RK2(pIn, pOut, t, dt, &parameters);
        tmp = pIn;
        pIn = pOut;
        pOut = tmp;
       
    }
    printf("t[%d] = %.16lf x = %.16lf v = %.16lf\n", i, (double) (t0+dt*(numberofSteps)), pOut->x, pOut->v);
   
   
   
    /*fp = fopen("data.dat", "w+");
    if ( ( fp=fopen("data.dat", "w+") ) == NULL) {
        printf("Couldn't open the file...Terminating\n");
        exit(1);
    }
    for (i = 0; i <= numberofSteps; ++i) {
       
        Point *tmp;
        double t = t0 + dt*i;
        fprintf(fp, "%.16lf %.16lf %.16lf\n", t, pIn->x, pIn->v);
       
        if (t == 80.) break;
       
        RK2(pIn, pOut, t, dt, &parameters);
        tmp = pIn;
        pIn = pOut;
        pOut = tmp;
       
    }
    fclose(fp);*/

   
   

  return 0;
}



/**************PROBLEM FUNCTIONS***************/

void drift(Point *pIn, Point *pOut, double t, Params *parameters) {
   
    //printf("n = %d\n", parameters->n);
    pOut->x = pIn->v;
    pOut->v = - pow(pIn->x, parameters->n) - ( 2. / t ) * pIn->v;
}



void RK2(Point *pIn, Point *pOut, double t, double dt, Params *parameters) {
   
    /* some temporary structs for help */
    double new_dt = 0.;
    Point tmp, tmp2;
    Point K1, K1dt;
    Point K2, K2dt;
    Point result;
   
    /* RungeKutta2 Implementation */
   
    //printP(pIn); 
    drift (pIn, &K1, t, parameters);
    //printP(&K1);
    HPMul(dt, &K1);
    pCpy(&K1, &K1dt);
    //printP(&K1dt);


    new_dt = dt * 0.5 ;
    pCpy(&K1, &tmp);
    //printP(&tmp);
    HPMul(0.5, &tmp);
    //printP(&tmp);
    drift(&tmp, &K2, t+new_dt , parameters);
    //printP(&K2);
    pCpy(&K2, &K2dt);
    HPMul(dt, &K2dt);
    //printP(&K2dt);
   
   
    pCpy(&K2dt, &result);
    pSum(&result, pIn);
    //printP(&result);
    pCpy(&result, &tmp2);
    //printP(&tmp2);
    pCpy(&tmp2, pOut);
    //printP(pOut);
   
   
}




/*********STRUCT FUNCTIONS (Tested beforehand)********/


void pSum(Point *res, Point *sum) {
   
    res->x = res->x + sum->x;
    res->v = res->v + sum->v;
}

void pMul(double factor, Point *res) {
   
    res->x = factor * res->x;
    res->v = factor * res->v;
}

void pFma(double factor, Point *res, Point *p1, Point *p2) { // fma(factor, p1, p2)  --> res = (factor * p1) + p2 at higher precision
    res->x = fma(factor, p1->x, p2->x);
    res->v = fma(factor, p1->v, p2->v);
}

void HPMul(double factor, Point *res) {
   
    Point add;
    add.x = 0;
    add.v = 0;
    Point *Null;
    Null = &add;
   
    res->x = fma(factor, res->x, Null->x);
    res->v = fma(factor, res->v, Null->v);
   
}

void pCpy(Point *pIn, Point *pOut) {
    pOut->x = pIn->x;
    pOut->v = pIn->v;
}

void pSet(double val1, double val2, Point *pIn) {
   
    if ( val2 == 0. ) {
        pIn->x = val1;
        pIn->v = val1;
    }
   
    else {
        pIn->x = val1;
        pIn->v = val2;
    }
   
}

void printP(Point *pIn) {
    printf("x  = %.8lf \t v = %.8lf\n", pIn->x, pIn->v) ;
}

Questo il plot per $n = 0 \to 4$

[SteezyMenchi]

Qui il testo dell'esercizio

[SteezyMenchi]

E qui il il metodo che ho implementato, preso da Numerical recipes in C (potrei aver sbagliato il nome ma tanto è conosciutissimo da quanto dice il professore)

[feddy]

Ciao,

non ho guardato il codice, ma sicuramente c'è un errore nell'algoritmo stando al plot che hai riportato. In questi casi conviene accertarsene con una soluzione analitica e verificare che sia riprodotta. Ad esempio, se integri con RK2 $y'=y, y(0)=1$, trovi come soluzione l'esponenziale $e^x$?
Nel tuo caso puoi prendere il caso vettoriale, magari ti è più comodo, cioè dato il vettore soluzione $y=(y_1,y_2)$, il tuo sistema è:
$$y_1' = y_1, y_2'=y_2$$ e dato iniziale il vettore $$y_0=(y_1(0),y_2(0))=(1,1)$$ che corrisponde al sistema $y'=Iy$.

Ricordati che in questi casi, per essere certi di aver implementato correttamente un metodo, basta calcolare l'ordine di convergenza e verificare che sia quello atteso (in questo caso 2).

[apatriarca]

Probabilmente al tuo professore non interessa, ma manca la verifica degli argomenti al programma. Dovresti almeno verificare che il numero di argomenti sia corretto.

Capisco che tu sia in qualche modo legato alla vecchia notazione, ma quando si fanno esercizi di questo tipo è a mio parere utile cercare di usare una notazione il più possibile simile al testo dell'esercizio. È soprattutto utile nella parte di visualizzazione degli output o richiesta di input. Per cui tu stai visualizzando il seguente testo:

Codice:

N = 0, t0 = 0.010000, x0 = 1.000000, v0 = 0.000000, dt = 0.010000, tmax = 100.000000
t[10000] = 100.0000000000000142 x = 1.5005478528675606 v = -100.0497285075012428


Tuttavia nel linguaggio del tuo esercizio t0, dt e tmax dovrebbero essere r0, dr e rmax. Similmente, x0 e v0 dovrebbero essere P(r0) e D(r0).

La variabile numberofSteps dovrebbe essere intera. Inoltre l'espressione usata per ottenere il risultato è errata. Un primo indizio sul fatto che sia errata è l'assenza di t0 dal calcolo. Se infatti passassi t0=99, dt=0.1 e tmax=100 avresti numberofSteps=1000 che è chiaramente sbagliato (sono necessari solo 10 step per arrivare da 99 a 100 con un passo da 0.1). Inoltre né il commento né l'espressione è chiara in quello che stai cercando di calcolare. Quel tipo di espressioni ha senso quando si ha a che fare con numeri interi, ma qui hai dei numeri in virgola mobile ed è più semplice usare la funzione corretta tra lceil, lfloor, ltrunc, lround... Suppongo tu voglia usare il valore intero più vicino per cui lround sarebbe quello corretto da usare.

Siccome il C supporta i numeri complessi avrei personalmente usato quelli invece di definire delle strutture per rappresentare dei vettori bidimensionali. In questo modo potevi scrivere espressioni usando la normale convenzione matematica.

Sinceramente farei a meno di cercare di generalizzare il codice più del necessario. In particolare l'uso della struttura Params sembra superflua. È certamente possibile scrivere una versione generica per RK ma è oltre le tue competenze e probabilmente sarebbe in generale più lenta di una versione specifica al tuo problema.

È inoltre possibile restituire strutture dalle funzioni e passarle come argomenti senza doverle passare per riferimento usando dei puntatori. In casi come questo in cui le strutture sono piccole è anzi forse anche desiderabile. Ci sono inoltre delle funzioni completamente superflue come quella per fare una copia o assegnare dei valori. Le strutture in C si possono infatti per esempio copiare semplicemente usando l'operatore di assegnamento:

Codice:

// pCpy(&tmp2, pOut) è uguale a
*pOut = tmp2;


Si possono inizializzare i valori di una struttura usando la seguente sintassi:

Codice:

Point zero = { 0.0, 0.0 };


e settare più valori in contemporanea come segue

Codice:

*pIn = (Point){ pIn->x + sum->x, pIn->v + sum->v };


Il C supporta inoltre nomi di funzioni più lunghi e credo che sia utile per incrementare la leggibilità del codice. Come piccolo esempio la funzione

Vedo che hai del codice per scrivere l'output su file invece che nella console. Credo sia utile sapere che è possibile farlo senza dover riscrivere il tuo programma con la seguente sintassi (operatore > seguito dal nome del file alla fine della riga di comando) in praticamente ogni sistema:

Codice:

./runge 3 0.01 1 0 0.01 100 > data.dat


Ovviamente questo lo puoi fare solo se i due output sono nello stesso formato.

L'esercizio ti chiede di stampare ogni tripletta (ri, Pi, Di) e non solo quella finale. C'è una ragione per cui hai commentato la riga nel ciclo per stamparle tutte? Nella consegna dell'esercizio devi certamente togliere il commento. C'è una ragione per cui stai usando 16 invece di 10 cifre decimali come richiesto dall'esercizio? Come puoi vedere le ultime cifre decimali mostrano principalmente errori numerici per cui sono praticamente inutili da conoscere.

Detto tutto questo non ti saprei ancora dire cosa ci sia esattamente di sbagliato nel codice perché la logica di Runge-Kutta appare molto offuscata e devo darci un occhiata più attenta.

[apatriarca]

Non è necessario usare nomi diversi per la struttura quando usi il typedef. Il seguente codice
è infatti perfettamente valido.

Codice:

struct Point {
    double x;
    double v;
};
typedef struct Point_t Point;


La funzione pSum, usata una singola volta, è equivalente al seguente codice:

Codice:

// pSum(&result, pIn);
result.x += pIn->x;
result.y += pIn->y;


Hai scritto più codice per definire la funzione di quello che va a sostituire.

Le funzioni pMul, pFma e pSet non sono state usate. Quindi possono essere rimosse.

La funzione printP è usata per debug. Il mio consiglio è quello di imparare ad usare un debugger, ma
abbiamo tutti fatto uso di funzioni di questo tipo e di printf per cui è forse la funzione che ha più senso tra quelle
che hai dichiarato.

La funzione HPMul è stata finalmente usata un po' di volte nel codice. Il codice è tuttavia abbastanza strano.
Inizi creando una struttura nulla, ne ottieni l'indirizzo e utilizzi quindi questo puntatore per accedere ai membri
della struttura originale che sappiamo essere uguali a zero. Di fatto il codice finale poteva essere scritto senza tutta
quella parte o l'uso di fma. Ho inserito alcuni commenti su possibili trasformazioni di alcune delle parti.

Codice:

// Point add;
// add.x = 0;
// add.v = 0;
// => Point add = {0, 0};
// Point *Null;
// Null = &add;
// => Point *Null = &add;
// res->x = fma(factor, res->x, Null->x);
// => res->x = fma(factor, res->x, 0);
// => res->x = factor * res->x;
res->x *= factor;
// res->v = fma(factor, res->v, Null->v);
res->y *= factor;


La funzione pCpy corrisponde ad una singola riga di codice che è più corta della chiamata alla funzione.
Quindi direi che non è necessaria.

A questo punto RK2 è diventata la seguente funzione. Qui il compilatore mi mostra un warning dicendo che K1dt è stata
settata ma non utilizzata. Non è un vero errore, ma mostra come il compilatore sia più capace di ragionare se non
facciamo uso di funzioni come quelle di cui hai fatto uso.

Codice:

double new_dt = 0.;
Point tmp, tmp2;
Point K1, K1dt;
Point K2, K2dt;
Point result;

drift(pIn, &K1, t, parameters);
// HPMul(dt, &K1);
K1.x *= dt;
K1.v *= dt;
// pCpy(&K1, &K1dt);
K1dt = K1;

new_dt = dt * 0.5 ;
// pCpy(&K1, &tmp);
tmp = K1;
// HPMul(0.5, &tmp);
tmp.x *= 0.5;
tmp.v *= 0.5;
drift(&tmp, &K2, t+new_dt , parameters);
// pCpy(&K2, &K2dt);
K2dt = K2;
// HPMul(dt, &K2dt);
K2dt.x *= dt;
K2dt.v *= dt;

// pCpy(&K2dt, &result);
result = K2dt;

// pSum(&result, pIn);
result.x += pIn->x;
result.v += pIn->v;

// pCpy(&result, &tmp2);
tmp2 = result;
// pCpy(&tmp2, pOut);
*pOut = tmp2;


È ora anche per noi più facile vedere un errore nella logica di Runge-Kutta. La variabile tmp è uguale
a 0.5 * K1 ma dovrebbe essere pIn + 0.5 * K1.

Una versione corretta del codice potrebbe forse essere qualcosa come segue:

Codice:

Point K1;
drift(pIn, &K1, t, parameters);

double half_dt = 0.5 * dt;
Point mid = {
    pIn->x + half_dt * K1.x,
    pIn->v + half_dt * K1.v
};

Point K2;
drift(&mid, &K2, t + half_dt, parameters);
pOut->x = pIn->x + dt * K2.x;
pOut->v = pIn->v + dt * K2.v;


o dopo aver cambiato le due funzioni per passare tutto per copia:

Codice:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

struct Point {
    double x;
    double v;
};
typedef struct Point Point;

Point drift(Point p, double t, double n);
Point RK2(Point p, double t, double dt, double n);

int main(int argc, char *argv[]) {
    if (argc != 7) {
        return 1;
    }

    int N = atoi(argv[1]);
    double t0 = atof (argv[2]);
    double x0 = atof(argv[3]);
    double v0 = atof(argv[4]);
    double dt = atof(argv[5]);
    double tmax = atof(argv[6]);
   
    int numberofSteps = (int) lround((tmax - t0) / dt);   // approximate to the nearest bigger integer

    Point p = { x0, v0 };
   
    printf("%.10f, %.10f, %.10f\n", t0, p.x, p.v);
    for (int i = 0; i < numberofSteps; ++i) {
        double t = t0 + dt*i;
        p = RK2(p, t, dt, N);
        printf("%.10f, %.10f, %.10f\n", t+dt, p.x, p.v);
    }
   
    return 0;
}

Point drift(Point p, double t, double n) {
    return (Point){
        p.v,
        - pow(p.x, n) - (2 / t) * p.v
    };
}

Point RK2(Point p, double t, double dt, double n) {
    Point K1 = drift(p, t, n);

    double half_dt = 0.5 * dt;
    Point mid = {
        p.x + half_dt*K1.x,
        p.v + half_dt*K1.v
    };

    Point K2 = drift(mid, t + half_dt, n);
    return (Point){
        p.x + dt*K2.x,
        p.v + dt*K2.v
    };
}


[apatriarca]

Alcuni problemi sono ovviamente dovuti alla tua familiarità con il linguaggio, ma credo che il problema principale sia che ti sei complicato la vita. Credo tu abbia cercato di usare degli strumenti di cui non hai la completa comprensione per principio piuttosto che per necessità. Hai insomma deciso che era necessario avere una funzione per calcolare le derivate, una per il calcolo di Runge-Kutte e delle strutture per poter passare gli argomenti a queste funzioni. E una volta che hai creato le strutture hai pensato a modi di usarle.. Vedi che nessuno di questi passaggi è effettivamente legato al tentare di risolvere il tuo problema che era di creare un programma che usasse RK2 per risolvere la tua equazione differenziale data. Da nessuna parte si è fatto riferimento a funzioni o strutture o altro. Il metodo più semplice era quindi quello di scrivere un programma che non facesse uso di funzioni o strutture o puntatori o ... Qualcosa come il seguente per esempio:

Codice:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int main(int argc, char *argv[argc+1]) {


    int n = atoi(argv[1]);
    double r0 = atof(argv[2]);
    double P0 = atof(argv[3]);
    double D0 = atof(argv[4]);
    double dr = atof(argv[5]);
    double rmax = atof(argv[6]);

    long steps = lround((rmax - r0) / dr);
    double half_dr = 0.5 * dr;

    double ri = r0;
    double Pi = P0;
    double Di = D0;

    printf("%.10f, %.10f, %.10f\n", ri, Pi, Di);

    for (long i = 0; i < steps; ++i) {       
        // Calcolo le derivate al punto ri
        double dPi = Di;
        double dDi = -pow(Pi, n) - (2/ri)*Di;

        // Calcolo il punto intermedio su cui calcolare le nuove derivate
        double r_mid = ri + half_dr;
        double P_mid = Pi + half_dr*dPi;
        double D_mid = Di + half_dr*dDi;

        // Calcolo le derivate al punto r_mid
        double dP_mid = D_mid;
        double dD_mid = -pow(P_mid, n) - (2/r_mid)*D_mid;

        // Aggiorno ri, Pi, Di
        ri += dr;
        Pi += dr*dP_mid;
        Di += dr*dD_mid;

        printf("%.10f, %.10f, %.10f\n", ri, Pi, Di);       
    }

    return 0;
}


Il codice è molto più corto e leggibile del tuo codice iniziale. Una volta che hai scritto qualcosa del genere puoi cercare di migliorarlo. Puoi per esempio osservare che il codice si ripete in modo abbastanza simile due volte. Questo è un buon segno che probabilmente vale la pena di usare una funzione e probabilmente ci sarà bisogno di strutture o altro. È comunque importante cercare di arrivare prima di tutto ad una soluzione e poi si generalizza o astrae. Se parti cercando di trovare una soluzione generale perdi di vista il tuo obiettivo e poi fai fatica a mettere insieme le cose.