// Einfuehrung in die Numerik, WS 2008
// Musterloesung Programmieraufgabe 1
// Gerd Ruecker
// Rundungsfehler beim Numerischen Differenzieren mit Differenzenquotienten 

#include <math.h>                 // Trigonometrische Funktionen
#include <stdio.h>                // Ausgabe in Datei
#include <stdlib.h>               // exit
#include <float.h>                // Bestimmung der Maschinengenauigkeit
#define FORMAT "%25.22f "         // Format der Float Ausgabe
#define STEPSIZE 0.08             // Logarithmische Schrittweite 
#define N_X 2                     // Anzahl der unabhaengigen Variablen
#define N_F 3                     // Anzahl der abhaengigen Variablen 

// Funktionen
double fcn1( double *x ){ return x[0]*x[0]*x[1] + 0.1 * exp(2* x[0] + x[1] ); }
double fcn2( double *x ){ return x[0] * x[1] + sin( x[0] * x[0] ); }
double fcn3( double *x ){ return sqrt( x[0] ) / x[1]; }

// Ableitungen 
double dfcn1( double *x, double *d ){
    return ( 2*x[0]*x[1] + 2*0.1*exp( 2*x[0]+x[1] ) ) * d[0] + ( x[0]*x[0] 
                       + 0.1 * exp( 2*x[0] + x[1] ) ) *d[1];        
}
double dfcn2( double *x, double *d ){
    return (x[1] + 2*x[0] * cos( x[0] * x[0] ) ) * d[0] + ( x[0] ) *  d[1] ;          
}
double dfcn3( double *x, double *d ) {
  return (1 / ( 2* sqrt( x[0] ) * x[1] ) ) * d[0]
                - ( sqrt( x[0] ) / ( x[1] * x[1] ) ) * d[1]; 
}

// Array von Funktionspointern aller Test Funktionen
double (*fcn[])( double *x ) = { fcn1, fcn2, fcn3 };
double (*dfcn[])( double *x, double *d ) = { dfcn1, dfcn2, dfcn3 };

// einseitiger und zweiseitige Differenzenquotientenberechnung
double DiffQuot( double *x, double (*fcn) ( double* ), double h, double *d, int mode ){
    int i;
		double f_arg1[N_X], f_arg2[N_X], one_side, two_side;
    for( i = 0; i < N_X; i++ )
       f_arg1[i] = x[i] + h*d[i], f_arg2[i] = x[i] - h*d[i];
    one_side = ( fcn( f_arg1 ) - fcn( x ) ) / h ;             // einseitiger DQ
    two_side = ( fcn( f_arg1 ) - fcn( f_arg2 ) ) / ( 2 * h ); // zweiseitiger DG
    return  mode == 0  ? one_side : two_side;
}

// Hauptprogramm
int main( int argc, char **argv ){   
   register int i;
	 FILE *fp;
   double x[N_X], d[N_X], f_d_val[N_F], f_d_val_c[N_F], error[N_F];
   register double h = 1.0, step = 1.0, eps = DBL_EPSILON;
   if( !( fp = fopen( argv[1], "w" ) ) ) exit( -1 );
   x[0] = 2.0, x[1] = 3.0, d[0] = 2.3, d[1] = 3.3;   
   do{    // Iteration mit Schrittweite h
     h = pow( 10, -step );   
     fprintf( fp, FORMAT, log10( h ) );
     for( i = 0; i < N_F; i++ ){
       f_d_val[i] = DiffQuot( x, fcn[i], h, d, 1 );
       f_d_val_c[i] = dfcn[i]( x, d );
       error[i] = fabs( f_d_val[i] - f_d_val_c[i] );
       fprintf( fp, FORMAT, log10( error[i] ));
     }
     fputc( '\n', fp );     
     step += STEPSIZE;
   }while( h >= eps );
   if( fclose( fp ) ) exit( -1 );
   return 0;
}