#include <math.h>
#include <stdio.h>
// Parametri della simulazione. Unità di misura SI
const double
h0= 2000, // Quota iniziale
g= 9.8, // Accelerazione di gravità
rho= 1, // Densità dell'aria
m= 80, // Massa
S= .7, // Area
Cr= 1, // Coeff. di resistenza
tolleranza= 1e-6, // Serve per controllare l'errore dell'integratore
dt_ini= 1, // Passo iniziale
k= .5 * rho * S * Cr,
k1= sqrt(g*m/k), // Questa è la velocità limite
k2= sqrt(g*k/m),
tc= acosh(exp(k2/k1 * h0)) / k2; // Tempo di caduta
static double emPOS= 0, emVEL= 0; // Errori massimi
// ============================================================
// === STATO ===
struct STS
{
double p, v;
};
STS operator*(const double mul, const STS &s)
{
STS r; r.p= mul * s.p; r.v= mul * s.v;
return r;
}
STS operator+(const STS &s1, const STS &s2)
{
STS r; r.p= s1.p + s2.p; r.v= s1.v + s2.v;
return r;
}
// ============================================================
int c_DER= 0; // Contatore chiamate a DER()
// Calcolo delle derivate
STS DER(double x, const STS &st)
{
// Forza aerodinamica
double Fa= k * st.v * st.v; // Fa = k * v^2
// Forza dovuta all'accelerazione di gravità
double Fg= m * g;
STS der; c_DER++;
der.p= st.v; // Derivata della posizione = velocità
der.v= (Fa - Fg) / m; // Derivata della velocità = accelerazione
return der; // Ritorna p' e v'
}
// ============================================================
////////////////////////////////////////////////////////////////////////////////
// //
// Description: //
// The Runge-Kutta-Fehlberg method is an adaptive procedure for approxi- //
// mating the solution of the differential equation y'(x) = f(x,y) with //
// initial condition y(x0) = c. This implementation evaluates f(x,y) six //
// times per step using embedded fourth order and fifth order Runge-Kutta //
// estimates to estimate the not only the solution but also the error. //
// The next step size is then calculated using the preassigned tolerance //
// and error estimate. //
// For step i+1, //
// y[i+1] = y[i] + h * ( 25 / 216 * k1 + 1408 / 2565 * k3 //
// + 2197 / 4104 * k4 - 1 / 5 * k5 ) //
// where //
// k1 = f( x[i],y[i] ), //
// k2 = f( x[i]+h/4, y[i] + h*k1/4 ), //
// k3 = f( x[i]+3h/8, y[i]+h*(3/32 k1 + 9/32 k2) ), //
// k4 = f( x[i]+12h/13, y[i]+h*(1932/2197 k1 - 7200/2197 k2 //
// + 7296/2197 k3) ), //
// k5 = f( x[i]+h, y[i]+h*(439/216 k1 - 8 k2 + 3680/513 k3 - 845/4104 k4))//
// k6 = f( x[i]+h/2, y[i]+h*(-8/27 k1 + 2 k2 - 3544/2565 k3 //
// + 1859/4104 k4 - 11/40 k5) ) //
// x[i+1] = x[i] + h. //
// //
// The error is estimated to be //
// err = h*( k1 / 360 - 128 k3 / 4275 - 2197 k4 / 75240 + k5 / 50 //
// + 2 k6 / 55 ) //
// The step size h is then scaled by the scale factor //
// scale = 0.8 * | epsilon * y[i] / [err * (xmax - x[0])] | ^ 1/4 //
// The scale factor is further constrained 0.125 < scale < 4.0. //
// The new step size is h := scale * h. //
////////////////////////////////////////////////////////////////////////////////
#define max(x,y) ( (x) < (y) ? (y) : (x) )
#define min(x,y) ( (x) < (y) ? (x) : (y) )
#define ATTEMPTS 12
#define MIN_SCALE_FACTOR 0.125
#define MAX_SCALE_FACTOR 4.0
static STS Runge_Kutta(STS *y, double x, double h);
////////////////////////////////////////////////////////////////////////////////
// int Embedded_Fehlberg_4_5( double (*f)(double, double), double y[], //
// double x, double h, double xmax, double *h_next, double tolerance ) //
// //
// Description: //
// This function solves the differential equation y'=f(x,y) with the //
// initial condition y(x) = y[0]. The value at xmax is returned in y[1]. //
// The function returns 0 if successful or -1 if it fails. //
// //
// Arguments: //
// double *f Pointer to the function which returns the slope at (x,y) of //
// integral curve of the differential equation y' = f(x,y) //
// which passes through the point (x0,y0) corresponding to the //
// initial condition y(x0) = y0. //
// double y[] On input y[0] is the initial value of y at x, on output //
// y[1] is the solution at xmax. //
// double x The initial value of x. //
// double h Initial step size. //
// double xmax The endpoint of x. //
// double *h_next A pointer to the estimated step size for successive //
// calls to Embedded_Fehlberg_4_5. //
// double tolerance The tolerance of y(xmax), i.e. a solution is sought //
// so that the relative error < tolerance. //
// //
// Return Values: //
// 0 The solution of y' = f(x,y) from x to xmax is stored y[1] and //
// h_next has the value to the next size to try. //
// -1 The solution of y' = f(x,y) from x to xmax failed. //
// -2 Failed because either xmax < x or the step size h <= 0. //
// //
////////////////////////////////////////////////////////////////////////////////
// //
int Embedded_Fehlberg_4_5(STS y[], double x, double h, double xmax, double *h_next, double tolerance )
{
double scale;
STS temp_y[2];
double err;
double yy;
int i;
int last_interval = 0;
// Verify that the step size is positive and that the upper endpoint //
// of integration is greater than the initial enpoint. //
if (xmax < x || h <= 0.0) return -2;
// If the upper endpoint of the independent variable agrees with the //
// initial value of the independent variable. Set the value of the //
// dependent variable and return success. //
*h_next = h;
y[1] = y[0];
if (xmax == x) return 0;
// Insure that the step size h is not larger than the length of the //
// integration interval. //
h = min(h, xmax - x);
// Redefine the error tolerance to an error tolerance per unit //
// length of the integration interval. //
// tolerance /= (xmax - x);
// Integrate the diff eq y'=f(x,y) from x=x to x=xmax trying to //
// maintain an error less than tolerance * (xmax-x) using an //
// initial step size of h and initial value: y = y[0] //
temp_y[0] = y[0];
temp_y[1] = y[0];
double x_old= 0; int n_stp= 0;
while ( x < xmax ) {
// Valori ESATTI della velocità e della quota
double v= -k1 * tanh(x * k2);
double p= h0 - m/k*log(cosh(x * k2));
double ePOS= temp_y[1].p - p; emPOS= fmax(emPOS, fabs(ePOS));
double eVEL= temp_y[1].v - v; emVEL= fmax(emVEL, fabs(eVEL));
printf("%3d) t= %7.3f", ++n_stp, x);
printf(" h= %8.3f (%8.1e)", temp_y[1].p, ePOS);
printf(" v= %8.3f (%8.1e)", temp_y[1].v, eVEL);
STS der= DER(0, temp_y[1]); c_DER--; // Questa chiamata NON è richiesta dall'integratore
printf(" a= %9f", der.v);
printf(" dt= %.3f\n", x - x_old); x_old= x;
if(temp_y[1].p < 0) break;
scale = 1.0;
for (i = 0; i < ATTEMPTS; i++) {
err = fabs( Runge_Kutta(temp_y, x, h).p ); // Errore su POS
if (err == 0.0) { scale = MAX_SCALE_FACTOR; break; }
yy = (temp_y[0].p == 0.0) ? tolerance : fabs(temp_y[0].p);
scale = 0.8 * sqrt( sqrt ( tolerance * yy / err ) );
scale = min( max(scale,MIN_SCALE_FACTOR), MAX_SCALE_FACTOR);
if ( err < ( tolerance * yy ) ) break;
h *= scale; //printf("h scalato di %g\n", scale);
if ( x + h > xmax ) h = xmax - x;
else if ( x + h + 0.5 * h > xmax ) h = 0.5 * h;
}
if ( i >= ATTEMPTS ) { *h_next = h * scale; return -1; };
temp_y[0] = temp_y[1];
x += h;
h *= scale;
*h_next = h;
if ( last_interval ) break;
if ( x + h > xmax ) { last_interval = 1; h = xmax - x; }
else if ( x + h + 0.5 * h > xmax ) h = 0.5 * h;
}
y[1] = temp_y[1];
printf("\nErrori massimi\n");
printf(" POS: %g m, VEL = %g m/s\n", emPOS, emVEL);
printf("\nPasso medio= %.3f s, chiamate a DER: %d\n", x / n_stp, c_DER);
return 0;
}
////////////////////////////////////////////////////////////////////////////////
// static double Runge_Kutta(double (*f)(double,double), double *y, //
// double x0, double h) //
// //
// Description: //
// This routine uses Fehlberg's embedded 4th and 5th order methods to //
// approximate the solution of the differential equation y'=f(x,y) with //
// the initial condition y = y[0] at x = x0. The value at x + h is //
// returned in y[1]. The function returns err / h ( the absolute error //
// per step size ). //
// //
// Arguments: //
// double *f Pointer to the function which returns the slope at (x,y) of //
// integral curve of the differential equation y' = f(x,y) //
// which passes through the point (x0,y[0]). //
// double y[] On input y[0] is the initial value of y at x, on output //
// y[1] is the solution at x + h. //
// double x Initial value of x. //
// double h Step size //
// //
// Return Values: //
// This routine returns the err / h. The solution of y(x) at x + h is //
// returned in y[1]. //
// //
////////////////////////////////////////////////////////////////////////////////
static const double a2 = 0.25;
static const double a3 = 0.375;
static const double a4 = 12.0 / 13.0;
static const double a6 = 0.5;
static const double b21 = 0.25;
static const double b31 = 3.0 / 32.0;
static const double b32 = 9.0 / 32.0;
static const double b41 = 1932.0 / 2197.0;
static const double b42 = -7200.0 / 2197.0;
static const double b43 = 7296.0 / 2197.0;
static const double b51 = 439.0 / 216.0;
static const double b52 = -8.0;
static const double b53 = 3680.0 / 513.0;
static const double b54 = -845.0 / 4104.0;
static const double b61 = -8.0 / 27.0;
static const double b62 = 2.0;
static const double b63 = -3544.0 / 2565.0;
static const double b64 = 1859.0 / 4104.0;
static const double b65 = -11.0 / 40.0;
static const double c1 = 25.0 / 216.0;
static const double c3 = 1408.0 / 2565.0;
static const double c4 = 2197.0 / 4104.0;
static const double c5 = -0.20;
static const double d1 = 1.0 / 360.0;
static const double d3 = -128.0 / 4275.0;
static const double d4 = -2197.0 / 75240.0;
static const double d5 = 0.02;
static const double d6 = 2.0 / 55.0;
// //
////////////////////////////////////////////////////////////////////////////////
// //
static STS Runge_Kutta(STS *y, double x0, double h)
{
STS k1, k2, k3, k4, k5, k6;
double h2 = a2 * h, h3 = a3 * h, h4 = a4 * h, h6 = a6 * h;
k1 = DER(x0, *y);
k2 = DER(x0+h2, *y + h * b21 * k1);
k3 = DER(x0+h3, *y + h * ( b31*k1 + b32*k2) );
k4 = DER(x0+h4, *y + h * ( b41*k1 + b42*k2 + b43*k3) );
k5 = DER(x0+h, *y + h * ( b51*k1 + b52*k2 + b53*k3 + b54*k4) );
k6 = DER(x0+h6, *y + h * ( b61*k1 + b62*k2 + b63*k3 + b64*k4 + b65*k5) );
*(y+1) = *y + h * (c1*k1 + c3*k3 + c4*k4 + c5*k5);
return d1*k1 + d3*k3 + d4*k4 + d5*k5 + d6*k6;
}
// ===============================================================================
int main(void)
{
STS y[2]; y[0].p= h0; y[0].v= 0; // Stato INIZIALE
double h_next;
int err= Embedded_Fehlberg_4_5(
y, // On input y[0] is the initial value of y at x, on output y[1] is the solution at xmax
0, // initial value of x
dt_ini, // Initial step size
99999, // xmax
&h_next, // A pointer to the estimated step size for successive calls to Embedded_Fehlberg_4_5
tolleranza // tolerance
);
printf("\nTempo di caduta: %.3f s, velocita' limite: %.3f m/s\n", tc, k1);
printf("\nCodice d'errore: %d\n", err);
return 0;
}
