#include <math.h>
#include <stdio.h>
// Usare "PROBLEMA 1" per verificare p(t) = 4t^2 + 3t + 2
// Usare "PROBLEMA 2" per verificare il corpo in caduta libera
#define PROBLEMA 1
// Usare "POSVEL 1" per verificare la derivata della posizione
// Usare "POSVEL 2" per verificare la derivata della velocità
#define POSVEL 1
// 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
t= 5, // Tempo da usare
k= .5 * rho * S * Cr,
k1= sqrt(g*m/k), // Questa è la velocità limite
k2= sqrt(g*k/m);
// ==========================================================
// Memorizza in fun il valore della funzione
// e in der il valore della sua derivata
void PVA(const double t, double &fun, double &der)
{
#if PROBLEMA == 1
double p= 4*t*t + 3*t + 2; // POS
double v= 8*t + 3; // VEL
double a= 8; // ACC
#elif PROBLEMA == 2
double p= m/k*log(cosh(t * k2)); // POS
double v= k1 * tanh(t * k2); // VEL
double a= (m * g - k * v * v) / m; // ACC
#else
#error PROBLEMA può valere 1 o 2
#endif
#if POSVEL == 1 // Derivata di POS
fun= p; der= v;
#elif POSVEL == 2 // Derivata di VEL
fun= v; der= a;
#else // ERRORE
#error POSVEL può valere 1 o 2
#endif
}
// Calcolo DERIVATA. "Numerical Recipes 3rd edition" Ch. 5.7 pag. 255
/*
Returns the derivative of a function func at a point x by Ridders’ method of polynomial extrapolation.
The value h is input as an estimated initial stepsize; it need not be small, but rather should be an increment in x
over which func changes substantially. An estimate of the error in the derivative is returned as err.
*/
double CalcDer(const double t, double hh, double &err)
{
const int ntab= 9; // Sets maximum size of tableau.
const double con= 1.4, con2= con*con; // Stepsize decreased by CON at each iteration.
const double safe= 2.0; // Return when error is SAFE worse than the best so far.
double ans, a[ntab][ntab];
hh *= con; // Questa moltiplicazione seguita dalla divisione "hh /= con" del ciclo for rende microsopicamente diversa hh originale e microscopicamente diverso l'errore stimato
err= HUGE_VAL;
for(int i= 0; i < ntab; i++) { // Successive columns in the Neville tableau will go to smaller stepsizes and higher orders of extrapolation.
hh /= con;
double fun, der, sp, sm;
PVA(t + hh, fun, der); sp= fun;
PVA(t - hh, fun, der); sm= fun;
a[0][i]= (sp - sm) / (2.0*hh);
if(i == 0) continue;
double fac= con2;
for(int j= 1; j <= i; j++) { // Compute extrapolations of various orders, requiring no new function evaluations.
a[j][i]= (a[j-1][i] * fac - a[j-1][i-1]) / (fac - 1.0);
fac= con2 * fac;
// The error strategy is to compare each new extrapolation to one order lower, both at the present stepsize and the previous one.
double errt= fmax(fabs(a[j][i] - a[j-1][i]), fabs(a[j][i] - a[j-1][i-1])); //printf("%d %d %-20.15g %-20.15g +/- %.15g\n", i, j, hh, a[j][i], errt);
if(errt <= err) { // If error is decreased, save the improved answer.
err= errt; ans= a[j][i];
}
}
if(fabs(a[i][i] - a[i-1][i-1]) >= safe * err) break; // If higher order is worse by a significant factor SAFE, then quit early.
}
return ans;
}
int main(void)
{
double hh= .1, err;
double fun, der_esa; PVA(t, fun, der_esa); printf("Derivata ESATTA = %.15g, funzione= %g\n", der_esa, fun);
// Rapporto incrementale
double ri, ri_old= HUGE_VAL, err_old= HUGE_VAL;
for(double DT= .5; DT > 1e-15; DT /= 2) {
double fm, fp, der; PVA(t - DT, fm, der); PVA(t + DT, fp, der); ri= (fp - fm) / (2*DT); //printf("%18.15g %g\n", ri, ri-ri_old);
if(fabs(ri - ri_old) > err_old) { ri= ri_old; break; }
err_old= fabs(ri - ri_old); ri_old= ri;
}
printf("\nRapp. incrementale= %.15g, errore= %g\n", ri, ri - der_esa);
double der_calc= CalcDer(t, hh, err);
printf("\nMetodo di Ridders = %.15g, incertezza= %g, errore= %g\n", der_calc, err, der_calc - der_esa);
return 0;
}
