Integration Methods (Euler Cromer,Verlet,Runge-Kutta, etc.)

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

#define EULER 1
#define EULER_CROMER 2
#define MID_POINT 3
#define HALF_PACE 4
#define VERLET 5
#define AUTO_VERLET 6
#define PRE_COR 7
#define RUNGE_KUTTA 8
#define POSITIVE 9
#define NEVER_MIND 10

struct pos_vel{
  double x;
  double v;
} xAndv0, xAndv;

double readdouble( double, int );

double energia( double, double, double, struct pos_vel, struct pos_vel, FILE* );

void eulero( double, double, double, double, struct pos_vel, struct pos_vel, FILE*, FILE* );

void eulero_cromer( double, double, double, double, struct pos_vel, struct pos_vel, FILE*, FILE* );

void punto_centrale( double, double, double, double, struct pos_vel, struct pos_vel, FILE*, FILE* );

void mezzo_passo( double, double, double, double, struct pos_vel, struct pos_vel, FILE*, FILE* );

void verlet( double, double, double, double, struct pos_vel, struct pos_vel, FILE*, FILE* );

void verlet_auto( double, double, double, double, struct pos_vel, struct pos_vel, FILE*, FILE* );

void pred_corr( double, double, double, double, struct pos_vel, struct pos_vel, FILE*, FILE* );

void runge_kutta( double, double, double, double, struct pos_vel, struct pos_vel, FILE*, FILE* );

double readdouble( double data, int pos ){
  scanf("%lf", &data);
  if( pos == POSITIVE ){
    while( data<=0 ){
      printf("Il valore inserito è errato! Reinseriscilo :\n");
      scanf("%lf", &data);
    }
  }
  printf("Il valore appena inserito è %lf. \n \n", data );
  return data;
}

double energia( double k, double m, double dt, struct pos_vel xAndv0, struct pos_vel xAndv, FILE *e){
  fprintf( e, "%lf                     %lf \n", dt, ((0.5*m*pow(xAndv.v,2.)+0.5*k*pow(xAndv.x,2.))-(0.5*m*pow(xAndv0.v,2.)+0.5*k*pow(xAndv0.x,2.)))/(0.5*m*pow(xAndv0.v,2.)+0.5*k*pow(xAndv0.x,2.)) );
  dt+=0.01;
  return dt;
}

void eulero( double k, double m, double dt, double T, struct pos_vel xAndv0, struct pos_vel xAndv, FILE *e, FILE *fp ){
  int j;
  double i, v;
  for( j=0; j<7; j++ ){
    xAndv=xAndv0;
    for(i=0; i<T; i+=dt){
      v=xAndv.v;
      xAndv.v-=(k/m)*xAndv.x*dt;
      xAndv.x+=v*dt;
      if(j==0){
    fprintf(fp, "%lf       %lf       %lf\n", i, xAndv.x, xAndv.v);
      }
    }
    dt=energia(k, m, dt, xAndv0, xAndv, e);
  }
}

void eulero_cromer( double k, double m, double dt, double T, struct pos_vel xAndv0, struct pos_vel xAndv, FILE *e, FILE *fp ){
  int j;
  double i;
  for( j=0; j<7; j++ ){
    xAndv=xAndv0;
    for(i=0; i<T; i+=dt){
      xAndv.v-=(k/m)*xAndv.x*dt;
      xAndv.x+=xAndv.v*dt;
      if(j==0){
    fprintf(fp, "%lf       %lf       %lf\n", i, xAndv.x, xAndv.v);
      }
    }
    dt=energia(k, m, dt, xAndv0, xAndv, e);
  }
}

void punto_centrale( double k, double m, double dt, double T, struct pos_vel xAndv0, struct pos_vel xAndv, FILE *e, FILE *fp ){
  int j;
  double i;
  for( j=0; j<7; j++ ){
    xAndv=xAndv0;
    for(i=0; i<T; i+=dt){
      xAndv.v-=(k/m)*xAndv.x*dt;
      xAndv.x+=xAndv.v*dt-0.5*(k/m)*xAndv.x*pow(dt,2.);
      if(j==0){
    fprintf(fp, "%lf       %lf       %lf\n", i, xAndv.x, xAndv.v);
      }
    }
    dt=energia(k, m, dt, xAndv0, xAndv, e);
  }
}

void mezzo_passo( double k, double m, double dt, double T, struct pos_vel xAndv0, struct pos_vel xAndv, FILE *e, FILE *fp ){
  int j;
  double i;
  for( j=0; j<7; j++ ){
    xAndv.x=xAndv0.x;
    xAndv.v=xAndv0.v-0.5*(k/m)*xAndv.x*dt;
    for(i=0; i<T; i+=dt){
      xAndv.v-=(k/m)*xAndv.x*dt;
      xAndv.x+=xAndv.v*dt;
      if(j==0){
    fprintf(fp, "%lf       %lf       %lf\n", i, xAndv.x, xAndv.v);
      }
    }
    dt=energia(k, m, dt, xAndv0, xAndv, e);
  }
}

void verlet( double k, double m, double dt, double T, struct pos_vel xAndv0, struct pos_vel xAndv, FILE *e, FILE *fp ){
  int j;
  double i, x_1, x=0.;
  for( j=0; j<7; j++ ){
    x_1=xAndv0.x;
    x=x_1+xAndv0.v*dt;
    for(i=0; i<T; i+=dt){
      xAndv.x=2.*x-x_1-(k/m)*x*pow(dt,2.);
      xAndv.v=(xAndv.x-x_1)/(2*dt);
      x_1=x;
      x=xAndv.x;
      if(j==0){
    fprintf(fp, "%lf       %lf       %lf\n", i, xAndv.x, xAndv.v);
      }
    }
    dt=energia(k, m, dt, xAndv0, xAndv, e);
  }
}

void verlet_auto( double k, double m, double dt, double T, struct pos_vel xAndv0, struct pos_vel xAndv, FILE *e, FILE *fp ){
  int j;
  double i, x;
  for( j=0; j<7; j++ ){
    xAndv=xAndv0;
    for(i=0; i<T; i+=dt){
      x=xAndv.x+xAndv.v*dt-0.5*(k/m)*xAndv.x*pow(dt,2.);
      xAndv.v+=0.5*(-(k/m)*xAndv.x - (k/m)*x)*dt;
      xAndv.x=x;
      if(j==0){
    fprintf(fp, "%lf       %lf       %lf\n", i, xAndv.x, xAndv.v);
      }
    }
    dt=energia(k, m, dt, xAndv0, xAndv, e);
  }
}

void pred_corr( double k, double m, double dt, double T, struct pos_vel xAndv0, struct pos_vel xAndv, FILE *e, FILE *fp ){
  int j;
  double i, xp, x, x_1, vn;
  for( j=0; j<7; j++ ){
    x_1=xAndv0.x;
    x=x_1+xAndv0.v*dt;
    vn-=(k/m)*xAndv0.x*dt;
    for(i=0; i<T; i+=dt){
      xp=x_1+2*xAndv.v*dt;
      xAndv.v=vn+((-(k/m)*xp-(k/m)*x)/2)*dt;
      xAndv.x=x+((xAndv.v+vn)/2)*dt;
      if(j==0){
    fprintf(fp, "%lf       %lf       %lf\n", i, xAndv.x, xAndv.v);
      }
      x_1=x;
      x=xAndv.x;
      vn=xAndv.v;
    }
    dt=energia(k, m, dt, xAndv0, xAndv, e);
  }
}

void runge_kutta( double k, double m, double dt, double T, struct pos_vel xAndv0, struct pos_vel xAndv, FILE *e, FILE *fp ){
  int j;
  double i, xp, vp;
  for( j=0; j<7; j++ ){
    xAndv=xAndv0;
    for(i=0; i<T; i+=dt){
      xp=xAndv.v*dt;
      vp=-(k/m)*xAndv.x*dt;
      xAndv.x+=(xAndv.v+0.5*vp)*dt;
      xAndv.v+=-(k/m)*(xAndv.x+0.5*xp)*dt;
      if(j==0){
    fprintf(fp, "%lf       %lf       %lf\n", i, xAndv.x, xAndv.v);
      }
    }
    dt=energia(k, m, dt, xAndv0, xAndv, e);
  }
}

main(){
  double k, m, dt, T;
  int i, algorithm;
  FILE *e, *fp;
  if((e=fopen("energia.dat", "w"))==NULL){
    printf("Errore nell'apertura del file energia.dat.\n");
    exit(EXIT_FAILURE);
  }
  if((fp=fopen("moto.dat", "w"))==NULL){
    printf("Errore nell'apertura del file energia.dat.\n");
    exit(EXIT_FAILURE);
  }
  fprintf(fp, "#Tempo         Spazio         Velocità\n");
  fprintf(e, "#Distanza temporale         Energia\n");
  printf("Il programma simula il moto di un oscillatore armonico facendo scegliere all'utente il metodo di integrazione più opportuno./n");
  printf("Il problema sta nel studiare come si muove un corpo di massa m attaccato ad una molla di costante elastica k che, al tempo t=0s, ha posizione x0 e velocità v0.\n");
  printf("Per iniziare quindi, ho bisogno che tu inserisca alcune quantità.\n");
  printf("Inserisci un valore di k positivo.\nk = ");
  k=readdouble( k, POSITIVE );
  printf("Inserisci un valore di m positivo.\nm = ");
  m=readdouble( m, POSITIVE );
  printf("Inserisci la posizione inziale x0.\nx0 = ");
  xAndv0.x=readdouble( xAndv0.x, NEVER_MIND );
  printf("Inserisci la velocità iniziale v0.\nv0 = ");
  xAndv0.v=readdouble( xAndv0.v, NEVER_MIND );
  printf("Inserisci il passo dt positivo con cui vuoi andare avanti nel tempo.\ndt = ");
  dt=readdouble( dt, POSITIVE );
  printf("Inserisci il tempo T positivo in cui l'analisi si deve fermare.\nT = ");
  T=readdouble( T, POSITIVE );
  printf("Ora che hai immesso i dati, scegli l'algoritmo di integrazione, immettendo il numero che lo affianca.\n");
  printf("1. Eulero\n2. Eulero-Cromer\n3. Metodo del punto centrale\n4. Metodo del mezzo passo\n5. Verlet\n6. Verlet autosufficiente\n7. Metodo di predizione e correzione\n8. Runge-Kutta\nLa tua scelta è : ");
  scanf("%d", &algorithm );
  while( algorithm<1 || algorithm>8 ){
    printf("Il valore inserito è errato! Reinseriscilo :\n");
    scanf("%d", &algorithm);
  }
  if( algorithm == EULER ){
    eulero( k, m, dt, T, xAndv0, xAndv, e, fp );
  }
  if( algorithm == EULER_CROMER ){
    eulero_cromer( k, m, dt, T, xAndv0, xAndv, e, fp );
  }
  if( algorithm == MID_POINT ){
    punto_centrale( k, m, dt, T, xAndv0, xAndv, e, fp );
  }
  if( algorithm == HALF_PACE ){
    mezzo_passo( k, m, dt, T, xAndv0, xAndv, e, fp );
  }
  if( algorithm == VERLET ){
    verlet( k, m, dt, T, xAndv0, xAndv, e, fp );
  }
  if( algorithm == AUTO_VERLET ){
    verlet_auto( k, m, dt, T, xAndv0, xAndv, e, fp );
  }
  if( algorithm == PRE_COR ){
    pred_corr( k, m, dt, T, xAndv0, xAndv, e, fp );
  }
  if( algorithm == RUNGE_KUTTA ){
    runge_kutta( k, m, dt, T, xAndv0, xAndv, e, fp );
  }
}