/* Test various methods on special Prothero--Robinson equation 
	y' = - C * ( y - phi( t ) ) + phi'( t )
	y(0) = y_0
*/
#include <math.h>
#include <stdio.h>
#include <time.h>
#include <unistd.h>

/* Max and min */
#define MIN(a,b) ((a) <= (b) ? (a) : (b))
#define MAX(a,b) ((a) >= (b) ? (a) : (b))

double phi(
	double t
)
{
	return cos(t);
}

double phiprime(
	double t
)
{
	return -sin(t);
}

double stepExplicitEuler(
	double h, 
	double t, 
	double C,
	double u
)
{
	return u + h * ( - C * ( u - phi( t ) ) + phiprime( t ) );
}

double stepImplicitEuler(
	double h, 
	double t, 
	double C,
	double u
)
{
	return ( u + h * ( C * phi( t + h ) + phiprime( t + h ) ) ) / ( 1.0 + h * C );
}

/* Number of methods to be tested: forward Euler, backward Euler */
#define nM 2

/* Number of C coefficients to test (increasing by factor of 10)*/
#define nC 6

/* Number of runs to make (increasing number of time steps s) */
#define nS 6

int
main(
	int argc,
	char **argv
)
{
	/* Local variables; method m, coefficient number c, number of steps s */
	int i, m, c, s, steps[nS];
	/* Initial value y0, final time T */
	double y0 = 2.0, T = 10.0;
	double err[nM][nC][nS], cpu[nM][nC][nS], coeff[nC];

	/* CPU time */
	clock_t CPUbegin;
	double CPUtime;

	/* Ptr to fcn which advances method one step for special PR eqn */
	double ( *stepFP )( double h, double t, double c, double u );

	/* Set C coefficient to increase by factor of 10 per run */
	coeff[0] = 1;
	for( c = 1; c < nC; c++ ) coeff[c] = 10 * coeff[c - 1];

	/* Set number of steps to increase by factor of 10 per run */
	steps[0] = 10;
	for( s = 1; s < nS; s++ ) steps[s] = 10 * steps[s - 1];
	
	/* Loop through methods */
	for( m = 0; m < nM; m++ ) {

		/* Test explicit Euler */
		if( m == 0 ) stepFP = &stepExplicitEuler;
		/* Test implicit Euler */
		if( m == 1 ) stepFP = &stepImplicitEuler;

		/* Loop through coeff values */
		for( c = 0; c < nC; c++ ) {
			double C = coeff[c];

			/* Loop through number of time steps */
			for( s = 0; s < nS; s++ ) {
				/* Number and size of time steps */
				int j, S = steps[s];
				double u, v, h = T / ( double ) S;

				/* Time S steps of desired method, quitting if solution is exploding */
				CPUbegin = clock();

				/* Initial values */
				u = y0;

				for( j = 1; j < S; j++ ) {
					double t = ( double ) j * h;
					u = ( *stepFP )( h, t, C, u );
					if( fabs( u ) > 1.0e+32 ) {
						break;
					}
				}

		        cpu[m][c][s] = ( clock() - CPUbegin ) / ( ( double ) CLOCKS_PER_SEC );

				/* Measure errors, cutoff to avoid underflow */
				v = exp( - MIN( 100.0, C * T ) ) * ( y0 - phi( 0.0 ) ) + phi( T );
				err[m][c][s] = fabs( u - v );

			} /* End for( s = 0; s < nS; s++ ) { */
		} /* End for( c = 0; c < nC; c++ ) { */
	} /* End for( m = 0; m < nM; m++ ) { */

	/* Print tables of errors and CPU times */

	for( m = 0; m < nM; m++ ) {

		if( m == 0 ) fprintf( stdout, "\n%%Errors in explicit Euler\n" );
		if( m == 1 ) fprintf( stdout, "\n%%Errors in implicit Euler\n" );
		
		/* Table header: S, h, CPU for first run */
		for( s = 0; s < nS; s++ ) fprintf( stdout, " & $S = %d$ ", steps[s] );
		fprintf( stdout, "\\\\ \n" );
		fprintf( stdout, " $C$	" );
		for( s = 0; s < nS; s++ ) fprintf( stdout, " & $h = %g$ ", T / steps[s] );
		fprintf( stdout, "\\\\ \n" );
		for( s = 0; s < nS; s++ ) fprintf( stdout, " & $CPU = %g$ ", cpu[m][0][s] );
		fprintf( stdout, "\\\\ \n" );

		for( c = 0; c < nC; c++ ) {
			fprintf( stdout, " %g ", coeff[c] );

			for( s = 0; s < nS; s++ ) fprintf( stdout, " & %6.3g", err[m][c][s] );
			fprintf( stdout, "\\\\ \n" );

		} /* End for( c = 0; c < nC; c++ ) { */

	} /* End for( m = 0; m < nM; m++ ) { */
	
	return 0 ;
}