/*
 * @(#) rungekt.c - Differential equation solver by Runge-Kutta method 3.
 * (c) 1996 Ivan Maidanski <ivmai@chat.ru> http://ivmai.chat.ru
 * Freeware program source. All rights reserved.
 **
 * Language: ANSI C
 * Tested with: Watcom C compiler v10
 * Last modified: 1996-03-06 15:35:00 GMT+03:00
 */

/* equation: y'(x)=f(x,y)=phi(x)*psi(x)-g(x)*y(x) */

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

#define PI 3.1415926

const double x0=-1;
const double y0=8;
const double b=2*PI-1;
const double eps1=1e-6;
const double eps=1e-3;

const int maxviewlines=15;

int n;
double *tablexa,*tableya,*tableyf;

double g(double x)
{
 return (-sin(x+1));
}

double phi(double x)
{
 return (sin(x+1));
}

double psi(double x)
{
 return (cos(x+1)/6);
}

double f(double x,double y)
{
 return (phi(x)*psi(x)-g(x)*y);
}

double yfunc(double x)
{ /* result function */
 double p=1-cos(x+1);
 return (exp(p)*8+p/6);
}

double euler(double xk,double yk,double h)
{
 return (f(xk,yk)*h+yk);
}

double rungekutta(double xk,double yk,double h)
{ /* precise level 3 */
 double p=f(xk,yk)/3;
 return (f(h*2/3+xk,f(h/3+xk,yk+p)*h*2/3+yk)+p)*h*3/4+yk;
}

int analyzeprecise(double dy)
{ /* for Runge-Kutta method with precise level 3 */
 dy=fabs(dy)/7;
 if (dy>=eps)
  return (1);
  else if (dy*16<eps)
   return (-1);
   else return (0);
}

double rungekuttaauto(double xk,double yk,double *h)
{ /* with step length auto-select */
 double y1=rungekutta(xk,yk,*h);
 double y2=rungekutta(xk,yk,(*h)/2);
 while (analyzeprecise(y1-y2)>0)
 {
  (*h)/=2;
  y1=y2;
  y2=rungekutta(xk,yk,*h/2);
 }
 while (analyzeprecise(y1-y2)<0 && (*h)*2<=b-xk)
 {
  (*h)*=2;
  y2=y1;
  y1=rungekutta(xk,yk,*h);
 }
 return ((analyzeprecise(y1-y2)<=0) ? y1 : y2);
}

void allocatetable(double **table)
{
 if (!(*table=(double *)realloc((void *)*table,
                                   n*sizeof(double))))
 {
  fprintf(stderr," No memory!\n");
  exit(1);
 }
}

void calcfixedtable(void)
{
 double h=(b-x0)/n;
 double x=x0;
 int i;
 tableyf=0;
 allocatetable(&tableyf);
 tableyf[0]=rungekutta(x0,y0,h);
 for (i=1;i<n;i++)
  tableyf[i]=rungekutta(x+=h,tableyf[i-1],h);
}

void calcautotables(void)
{
 double h=(b-x0)/2;
 n=1;
 tablexa=tableya=0;
 allocatetable(&tablexa);
 allocatetable(&tableya);
 tableya[0]=rungekuttaauto(x0,y0,&h);
 tablexa[0]=x0+h;
 while (tablexa[n-1]<b-eps1)
 {
  n++;
  allocatetable(&tablexa);
  allocatetable(&tableya);
  tableya[n-1]=rungekuttaauto(tablexa[n-2],tableya[n-2],&h);
  tablexa[n-1]=tablexa[n-2]+h;
 }
}

int getnexttablexapos(double x,int k)
{
 for (k=k;k<n;k++)
  if (tablexa[k]>=x)
   return (k);
 return (-1);
}

double approximatebyline(double x,
                         double x1,double y1,double x2,double y2)
{
 return ((y2-y1)*(x-x1)/(x2-x1)+y1);
}

double gettableyawithapprox(double x,int *k)
{
 *k=getnexttablexapos(x,*k);
 if (tablexa[*k]==x)
  return (tableya[*k]);
  else return ((*k==0) ?
               approximatebyline(x,x0,y0,tablexa[0],tableya[0]) :
               approximatebyline(x,tablexa[*k-1],tableya[*k-1],
                                 tablexa[*k],tableya[*k]));
}

void printresulttable(void)
{
 int i,k=0;
 int m= (n>maxviewlines)? 1 : n/maxviewlines;
 double h=(b-x0)/n;
 double x;
 printf("....N.....-X-......-Yr-.....-Ya-.....-Yf-.\n");
 printf(" %5ld %8.3f %8.3f %8.3f %8.3f\n",0,x0,y0,y0,y0);
 for (i=0;i<n;i+=m)
  printf(" %5ld %8.3f %8.3f %8.3f %8.3f\n",i+1,x=h*i+x0,
         yfunc(x),gettableyawithapprox(x,&k),tableyf[i]);
 printf("\n");
}

/*
double getdeltamax(double *table)
{
 int i;
 double max=fabs(y0);
 for (i=0;i<n;i++)
  if (fabs(table[i])>max)
   max=fabs(table[i]);
 return (max);
}

double getdeltaapproxmax(double *table)
{
 int i;
 double val,max=0;
 double h=(b-x0)/n;
 double x=x0;
 for (i=0;i<n;i++)
  if ((val=fabs(yfunc(x+=h)-table[i]))>max)
   max=val;
 return (max);
}

void drawaxes(int dx,int dy)
{
 int xmax=getmaxx-1,ymax=getmaxy-1;
 setcolor(white);
 line(1,ymax/2,xmax,ymax/2);
 line(xmax/2,1,xmax/2,ymax);
 line(xmax,ymax/2,xmax-5,ymax/2-3);
 line(xmax,ymax/2,xmax-5,ymax/2+3);
 line(xmax/2,1,xmax/2-3,6);
 line(xmax/2,1,xmax/2+3,6);
 line(xmax/2+dx,ymax/2-3,xmax/2+dx,ymax/2+3);
 line(xmax/2-3,ymax/2-dy,xmax/2+3,ymax/2-dy);
}
*/

void main(void)
{
 printf("\n\n Differential Equation Solver\n "
        " by Runge-Kutta method (precise level 3)\n"
        " Equation:\n "
        " y'(x)=f(x,y)=phi(x)*psi(x)-g(x)*y(x),\n "
        " where g(x)=-sin(x+1), phi(x)=sin(x+1), psi(x)=cos(x+1)/6\n"
        " x0=%.3f, y0=%.3f,\n"
        " b(xn)=%.3f, eps1=%g\n\n",x0,y0,b,eps1);
 printf("Solving with auto-steps...\n");
 calcautotables();
 printf("Solving with %d fixed steps...\n",n);
 calcfixedtable();
 printf("Result table:\n");
 printresulttable();
 free((void *)tablexa);
 free((void *)tableya);
 free((void *)tableyf);
}