/*
 *  ODEの爆発時刻および爆発レートを求めるプログラム
 *  本プログラムでは 
 *    y'=y^2 ,  y(0)=1
 *  を解く．解は
 *    y(t)=1/(1-t)
 *  となり，t=1で爆発レート1で爆発する．
 *
 *  関数func と 初期条件さえ書き換えれば，他の問題にも使用可能
 *  
 *  コンパイルは以下のように行う
 *  % gcc blowup.c dopri5.c -lm
 *
 *  cf 廣田千明，小澤一文，常微分方程式系の解の爆発時刻および爆発レートの推定法
 *     -偏微分方程式の爆発問題への応用-，日本応用数理学会論文誌，14 (2004), 
 *     13-38.
 *
 *  2004.09.30 作成
 *  2004.10.25 修正  soutをグローバル変数に
 */

#include <math.h>
#include <stdio.h>
#include <string.h>
#include "dopri5.h"


#define  ndgl      2 /* 方程式の次元+1 */
#define  nrdens    2 /* dense output する成分の数 */
#define  licont    nrdens


#define b_time 1.0 /* blow-up time(誤差を計算するためにのみ利用) */
#define b_rate 1.0 /* blow-up rate(誤差を計算するためにのみ利用) */


/* 数列Sl = S0 * gamma^l, l=0,1,2,...,l_max  */
#define S0 16
#define gamma 2    
#define l_max 10 


#define n_stage 3 /* Aitken \Delta^2法の適用回数 */
double tl[n_stage+1][l_max+1];
int l_count=0;

unsigned int steps[l_max+1]; /* step数を格納する入れ物 */


/* 微分方程式の右辺のf */
void func (unsigned n, double t, double *y, double *f)
{
	int i;
	f[0] = y[0]*y[0];
}

/* 弧長変換を施したf(書き換え不要) */
void trns_func (unsigned n, double s, double *y, double *f)
{
	int i;
	double t,norm=0.0;

	t=y[n-1];

	func(n-1,t,y,f);
	f[n-1] = 1.0;

	for (i=0;i<n;i++)
		norm += f[i]*f[i];
	norm = sqrt(norm);

	for (i=0;i<n;i++)
		f[i] = f[i]/norm;
}

double sout;

/* 数値解を出力する関数(dopri5内で使用される) */
void solout (long nr, double sold, double s, double* y, unsigned n, int* irtrn)
{
  double y1;

  if (nr == 1)
  {
		printf ("      s=%9.3e t=%9.3e y=%9.3e step=%d\n", s, y[n-1], y[0], nr-1);
    sout = S0;
		steps[l_count]=nr-1;
  }
  else
    while (s >= sout)
			{
				printf ("l=%3d s=%9.3e t=%9.3e y=%9.3e step=%d\n",l_count, sout, contd5(n-1,sout), contd5(0,sout), nr-1);
				tl[0][l_count]=contd5(n-1,sout);
				steps[l_count]=nr-1;
				l_count++;
				sout = sout*gamma;
			}
	
}

int main (void)
{
  double   y[ndgl];
  unsigned icont[licont], i, l, stage;
  int      res, iout, itoler;
  double   s, t, s_end, atoler, rtoler;
	char format0[16],format1[256];
	double estimated_T,rate[l_max+1];

  iout = 2; /* dense output を行う */
  s = 0.0; /* 弧長 */
	t = 0.0; /* 時刻 */

  y[0] = 1.0; /* 初期条件 */ 

	y[ndgl-1] = t; 
  s_end = S0*pow(gamma,l_max);

	/* 許容誤差の指定 */
  itoler = 0;
  rtoler = 1.0E-15;
  atoler = rtoler;

	/* dense output を用いる成分(書き換え不要) */
	for(i=0;i<ndgl;i++)
		icont[i] = i;

  res = dopri5 (ndgl, trns_func, s, y, s_end, &rtoler, &atoler, itoler,
								solout, iout, stdout, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
								0, 0, 0, nrdens, NULL, licont);

	printf ("rtol=%9.3e   fcn=%li   step=%li   accpt=%li   rejct=%li\r\n",
					rtoler, nfcnRead(), nstepRead(), naccptRead(), nrejctRead());
		

	/* Aitken加速 */
	for (stage=1;stage<=n_stage;stage++)
		for (l=2*stage;l<=l_max;l++)
			tl[stage][l]=tl[stage-1][l]-pow(tl[stage-1][l-1]-tl[stage-1][l],2)
				/(tl[stage-1][l]-2.0*tl[stage-1][l-1]+tl[stage-1][l-2]);


	/* 爆発レートの計算 */
	estimated_T=tl[n_stage-1][l_max];
	for(l=1;l<=l_max;l++){
		rate[l]=-log(gamma)/log(fabs((tl[0][l]-estimated_T)
																 /(tl[0][l-1]-estimated_T)));
	}


	/* 結果の出力 */
	strcpy(format1,"  l  steps     rate   ");
	for(stage=0;stage<=n_stage;stage++){
		sprintf(format0,"     t%d    ",stage);
		strcat(format1,format0);
	}
	strcat(format1,"\n");

	/* 爆発時刻の出力 */
	puts("\n");
	printf("blow-up time\n");
	printf(format1);
	printf("%3d %6d             ",0,steps[0]);
	printf("%+10.3e \n",tl[0][0]);
	for(l=1;l<=l_max;l++) {
		printf("%3d %6d %+10.3e  ",l,steps[l],rate[l]);
		for(stage=0;stage<= ((n_stage<(l/2)) ? n_stage : (l/2));stage++)
			printf("%+10.3e ",tl[stage][l]);
		printf("\n");
	}
	
	/* 誤差の出力 */
	puts("\n");
	printf("error\n");
	printf(format1);
	printf("%3d %6d             ",0,steps[0]);
	printf("%+10.3e \n",tl[0][0]-b_time);
	for(l=1;l<=l_max;l++) {
		printf("%3d %6d %+10.3e  ",l,steps[l],rate[l]-b_rate);
		for(stage=0;stage<= ((n_stage<(l/2)) ? n_stage : (l/2));stage++)
			printf("%+10.3e ",tl[stage][l]-b_time);
		printf("\n");
	}

 	return 0;

}