Compare Euler's, Improved Euler's, and Runge-Kutta methods of Estimating Differential Equation

compare_e_ie_r.c

/**
Comparison between methods (euler, imp euler, & runge-kutta)
@author Tyler Weaver

Compile:
gcc compare.c -o compare
Run:
./compare

*/

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

#define X 	0
#define Y 	1

double func_Dy(double x, double y); // y' = (y + x)^2
double func_y(double x); // y = tan(x) - x
double runge_kutta(double h, double target_x, 
	double initial_x, double initial_y,
	double (*func)(double,double), 
	double (*actual)(double),
	bool debug);
double improved_euler(double step_size, double target_x, 
	double initial_x, double initial_y,
	double (*func)(double,double), 
	double (*actual)(double),
	bool debug);
double euler(double step_size, double target_x, 
	double initial_x, double initial_y,
	double (*func)(double,double), 
	double (*actual)(double),
	bool debug);

int main()
{
	const double h = 0.1;
	const double initial_value[2] = {0,1};
	const double target_x = 1.5;
	double output;

	euler(h, target_x, initial_value[X], initial_value[Y],
		func_Dy, func_y, true);
	improved_euler(h, target_x, initial_value[X], initial_value[Y],
		func_Dy, func_y, true);
	runge_kutta(h, target_x, initial_value[X], initial_value[Y],
		func_Dy, func_y, true);

	return 0;
}

double func_Dy(double x, double y)
{
	double value = x*y + y;
	return value;
}

double func_y(double x)
{
	double y = exp((x*(x + 2))/2);
	return y;
}

double runge_kutta(double h, 
	double target_x, double x, double y,
	double (*func)(double,double),
	double (*actual)(double), bool debug)
{
	static bool first = true;

	double target_y; // the output value we are looking for
	double exact, error;
	double k[4];

	if(first && debug) // first time through, print the table header
	{
		printf("Runge-Kutta method\r\n");
		printf("x_n,   y_n,     exact,   error\r\n");
		printf("------|--------|--------|------\r\n");
		first = false;
	}

	if(x >= target_x) // we have reached the end
		return y;
	else
	{
		k[0] = h*func(x,y);						// k1
		k[1] = h*func(x + .5*h, y + .5*k[0]); 	// k2
		k[2] = h*func(x + .5*h, y + .5*k[1]);	// k3
		k[3] = h*func(x + h, y + k[2]);			// k4

		x = x + h;				// update x
		y = y + (k[0] + 2*k[1] + 2*k[2] + k[3])/6.0;

		exact = actual(x);
		error = fabsf(exact - y);

		if(debug)
			printf("%3.1f, %7.3f, %7.3f, %7.3f\r\n", x,y,exact,error); // output values for debuging
		target_y = runge_kutta(h,target_x,x,y,func,actual,debug); // do it again...
	}
	first = true; // reset the first flag
	return target_y; // return the target_y value
}

double improved_euler(double step_size, 
	double target_x, double x, double y,
	double (*func)(double,double),
	double (*actual)(double), bool debug)
{
	static bool first = true;

	double target_y; // the output value we are looking for
	double exact, error;
	double predictor;

	if(first && debug) // first time through, print the table header
	{
		printf("Improved Euler's method\r\n");
		printf("x_n,   y_n,     exact,   error\r\n");
		printf("------|--------|--------|------\r\n");
		first = false;
	}

	if(x >= target_x) // we have reached the end
		return y;
	else
	{
		predictor = y + (step_size*func(x, y)); // calculate the predictor
		// now the corrector
		y = y + 0.5*step_size*(func(x,y) + func(x+step_size,predictor)); 
		x = x + step_size;				// update x
		exact = actual(x);
		error = fabsf(exact - y);
		if(debug)
			printf("%3.1f, %7.3f, %7.3f, %7.3f\r\n", x,y,exact,error); // output values for debuging
		target_y = improved_euler(step_size,target_x,x,y,func,actual,debug); // do it again...
	}
	first = true; // reset the first flag
	return target_y; // return the target_y value
}

double euler(double step_size, 
	double target_x, double x, double y,
	double (*func)(double,double),
	double (*actual)(double), bool debug)
{
	static bool first = true;

	double target_y; // the output value we are looking for
	double exact, error;

	if(first && debug) // first time through, print the table header
	{
		printf("Euler's method\r\n");
		printf("x_n,   y_n,     exact,   error\r\n");
		printf("------|--------|--------|------\r\n");
		first = false;
	}

	if(x >= target_x) // we have reached the end
		return y;
	else
	{
		y = y + (step_size*func(x, y)); // calculate the new y value
		x = x + step_size;				// update x
		exact = actual(x);
		error = fabsf(exact - y);
		if(debug)
			printf("%3.1f, %7.3f, %7.3f, %7.3f\r\n", x,y,exact,error); // output values for debuging
		target_y = euler(step_size,target_x,x,y,func,actual,debug); // do it again...
	}
	first = true; // reset the first flag
	return target_y; // return the target_y value
}