/**************************************************************************************
 * Numerical Analysis Programming #1 : 4th order Runge-Kutta Method
 * ------------------------------------------------------------------------------------
 * PROBLEM : Solve x' = 2+(x-t-1)^2 with x(1)=2 on the interval [1,1.5625] using 
 *           Runger-Kutta Method. The step size is calculated by dividing the length of 
 *           the interval by the number of steps,n=72. Exact solution x(t)=1+t+tan(t-1) 
 *
 * INPUT   : Initial value, The number of steps, Interval  
 * OUTPUT  : Numerical Solution, Exact Solution, Error 
 *
 * VARIABLE DECLARATIONS 
 *     1. double h : step size
 *     2. double t : mesh points over the interval
 *     3. double x : numerical solution at t
 * ------------------------------------------------------------------------------------
 * Source code File: ~kkjin/text-340-rk4/rk4.c
 * Executable  File: ~kkjin/text-340-rk4/rk4
 **************************************************************************************/

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

 #define   LOWER_B   1.0       /* lowerbound of interval [1,1.5625] */
 #define   UPPER_B   1.5625    /* upperbound of interval [1,1.5625] */
 #define   STEPS     72        /* number of steps */
 #define   INIT      2.0       /* initial value x(1)=2 */

 double    F(double,double);         /* Function Prototype Declarations */
 double    Exact_Solution(double);

 main()
 {
        double h,x,t;                /* Variable Declarations */
	double f1,f2,f3,f4,error,exact;
	int    i; 

        t=LOWER_B; x=INIT; i=1;      /* Initialization */
	h=(UPPER_B-LOWER_B)/STEPS;
	
        printf(" ====================================================================\n");
        printf("   n          t              x(t)          Real Sol          Error\n");
        printf(" --------------------------------------------------------------------\n");

        while(t<UPPER_B){             
           f1 = h*F(t,x);                    /* 4th Order Runge-Kutta Formula */
           f2 = h*F(t+h/2.0,x+f1/2.0);
           f3 = h*F(t+h/2.0,x+f2/2.0);
           f4 = h*F(t+h,x+f3);
           
           x += (f1+2*(f2+f3)+f4)/6.0;
	   t  = LOWER_B + h*i;               /* update t */

           exact = Exact_Solution(t);        /* find exact solution at t */
           error = fabs(x-exact);            /* calculate error */
           printf("%4d %15.8lf %15.8lf %15.8f %15.8lf\n",i++,t,x,exact,error);
        }
 }




/* Function : double F(double,double)
 * -----------------------------------------
 *  Calculate f(t,x) and return the result. 
 */
 double F(double t,double x)
 {
        double value;
        value = 2+pow(x-t-1,2);

	return value; 
 }



 
/* Function : double Exact_Solution(double)
 * ----------------------------------------
 *  Calculate exact solution of given differential equation at t.
 */
 double Exact_Solution(double t)
 {
        double exact;
        exact = 1+t+tan(t-1);
	
	return exact;
 }