/*************************************************************************
*				PECE6.CPP                                *
*		METHODE DE PREDICTION CORRECTION D'ORDRE 6               *
*            pour l'integration d'une equation differentielle            *
*	      		   de la forme Y'=F(x,Y)                         *
*************************************************************************/
/*	 		AMELIORATIONS A APPORTER
 integration sur l'intervalle [0,1] uniquement, à modifier. De plus
la fonction F(x,Y) devra etre calculee a part, dans une fonction separee*/

#include<iostream.h>
#include<math.h>
#include<conio.h>
#include<stdio.h>
#include<stdlib.h>
main()
{
clrscr();
int i;
unsigned long n,nbpoint;
double Y0,Y,F[8],F0,AB[7][7],AM[7][7],AMK[7],h,somme,som,x;
//**********CREATION  DU  FICHIER  DE  RESULTATS : RESULT.DAT*************
FILE *res;
res=fopen("c:\\tc\\bin\\result.dat","w");

//**********initialisation et coefficients********************************
Y=1;

AB[0][0]=1;
AB[0][1]=(double)1.5;
AB[1][1]=(double)-.5;
AB[0][2]=(double)23/12;
AB[1][2]=(double)-4/3;
AB[2][2]=(double)5/12;
AB[0][3]=(double)55/24;
AB[1][3]=(double)-59/24;
AB[2][3]=(double)37/24;
AB[3][3]=(double)-3/8;
AB[0][4]=(double)1901/720;
AB[1][4]=(double)-1387/360;
AB[2][4]=(double)109/30;
AB[3][4]=(double)-637/360;
AB[4][4]=(double)251/720;
AB[0][5]=(double)4277/1440;
AB[1][5]=(double)-7923/1440;
AB[2][5]=(double)4991/720;
AB[3][5]=(double)-3649/720;
AB[4][5]=(double)959/480;
AB[5][5]=(double)-95/288;
AB[0][6]=(double)198721/60480;
AB[1][6]=(double)-18637/2520;
AB[2][6]=(double)235183/20160;
AB[3][6]=(double)-10754/945;
AB[4][6]=(double)135713/20160;
AB[5][6]=(double)-5603/2520;
AB[6][6]=(double)19087/60480;

AM[0][0]=.5;
AM[0][1]=(double)2/3;
AM[1][1]=(double)-1/12;
AM[0][2]=(double)19/24;
AM[1][2]=(double)-5/24;
AM[2][2]=(double)1/24;
AM[0][3]=(double)323/360;
AM[1][3]=(double)-11/30;
AM[2][3]=(double)53/360;
AM[3][3]=(double)-19/720;
AM[0][4]=(double)1427/1440;
AM[1][4]=(double)-133/240;
AM[2][4]=(double)241/720;
AM[3][4]=(double)-173/1440;
AM[4][4]=(double)3/160;
AM[0][5]=(double)2713/2520;
AM[1][5]=(double)-15487/20160;
AM[2][5]=(double)586/945;
AM[3][5]=(double)-6737/20160;
AM[4][5]=(double)263/2520;
AM[5][5]=(double)-863/60480;
AM[0][6]=(double)139849/120960;
AM[1][6]=(double)-121797/120960;
AM[2][6]=(double)123133/120960;
AM[3][6]=(double)-88547/120960;
AM[4][6]=(double)41499/120960;
AM[5][6]=(double)-11351/120960;
AM[6][6]=(double)275/24192;

AMK[0]=(double).5;
AMK[1]=(double)5/12;
AMK[2]=(double)9/24;
AMK[3]=(double)251/720;
AMK[4]=(double)95/288;
AMK[5]=(double)19087/60480;
AMK[6]=(double)36799/120960;


cout<<"\nnombre de points dans la discretisation :";
cin>>nbpoint;

h=(double)1/nbpoint;cout<<"\nh="<<h<<'\n';
//**********DEMARRAGE  DE  LA  METHODE************************************
//calcul des F[0] pour la premiere fois**************************************
		F[0]=Y;
//BOUCLES DE DEMARRAGE****************************************************
for(n=0;n<=6;n++)
	{
	x=(n+1)*h;
//prediction dans la phase de demarrage***********************************
		 somme=0;
		 for(i=0;i<=n;i++)
			 somme=somme+AB[i][n]*F[n-i];
		 Y0=Y+somme*h;
//evaluation dans la phase de demarrage***********************************
		F0=Y0+x;
//correction dans la phase de demarrage***********************************
		som=0;
		for(i=0;i<=n;i++)
			som=som+AM[i][n]*F[n-i];
		Y=Y+(som+AMK[n]*F0)*h;cout<<Y<<'\n';
//evaluation dans la phase de demarrage***********************************
		F[n+1]=Y+x;
	}
//reclassement dans la phase de demarrage*********************************
	for(i=0;i<=6;i++)
		F[i]=F[i+1];
//**********METHODE  D'ORDRE  6*******************************************
for(n=7;n<=nbpoint-1;n++)
	{
	x=(n+1)*h;
//prediction**************************************************************
		 somme=0;
		 for(i=0;i<=6;i++)
			 somme=somme+AB[i][6]*F[6-i];
		 Y0=Y+somme*h;
//evaluation**************************************************************
		F0=Y0+x;
//correction**************************************************************
		som=0;
		for(i=0;i<=6;i++)
			som=som+AM[i][6]*F[6-i];
		Y=Y+(som+AMK[6]*F0)*h;
//evaluation**************************************************************
		F[7]=Y+x;
//*****reclassement dans la methode d'ordre 6*****************************
		for(i=0;i<=6;i++)
			F[i]=F[i+1];
//resultats dans result.dat tous les 6 degres*****************************
	if((n+1)%10==0)
	  {
	  cout<<x<<'\t'<<Y<<'\t'<<2*exp(x)-x-1<<'\n';
	  fprintf(res,"%lf %lf %lf\n",x,Y,2*exp(x)-x-1);
	  }
	} //fin de la boucle sur n
getch();
return 0;
}