//ODE_Gear.cpp	積分剛性方程組吉爾法

#include <iostream>							//輸入輸出流頭文件
#include "OrdinaryDifferentialEguation.h"	//求解常微分方程(組)頭文件
using namespace std;						//名字空間


void main(void)
{	
	double a(0.0), b(1.0), h(0.01), hmax(0.1);
	double hmin[4]={0.0001,0.0001,0.00001,0.00001};
    double eps[4]={0.001,0.0001,0.0001,0.00001};
	matrix<double> z(3,30);
	valarray<double> y(3), t(30);
	
  	cout.setf(ios::fixed);	//輸出數據為定點法
	cout.precision(6);		//精度6位

	for(int k=0; k<4; k++)
	{
		y[0]=1.0;
		y[1]=0.0; 
		y[2]=-1.0;

		int iReValue = ODE_Gear(a,b,hmin[k],hmax,h,eps[k],y,30,t,z);
		cout << "h = " << h << endl;
		cout << "hmin = " << hmin[k] << endl;
		cout << "hmax = " << hmax << endl;
		cout << "eps = " << eps[k] << endl << endl;
		cout << "iReValue = " << iReValue << endl << endl;

		for(int i=0; i<30; i++)
		{
			cout << "t = " << t[i] << endl;
			for(int j=0; j<3; j++)
				cout << "    y(" << j << ") = " << z(j, i) << "\t";
			cout << endl;
		}
		cout << endl;
	}
}

//計算微分方程組中各方程右端函數值
template <class _Ty>
void FunctionValueG(_Ty t, valarray<_Ty>& y, valarray<_Ty>& d)
{
	int n = y.size();	//微分方程組中方程的個數，也是未知函數的個數

    d[0]=-21.0*y[0]+19.0*y[1]-20.0*y[2];
    d[1]=19.0*y[0]-21.0*y[1]+20.0*y[2];
    d[2]=40.0*y[0]-40.0*y[1]-40.0*y[2];

}

//計算雅可比矩陣的函數
template <class _Ty>
void JacobiMatrix(_Ty t, valarray<_Ty>& y, matrix<_Ty>& p)
{ 
	  int n = y.size();	//微分方程組中方程的個數，也是未知函數的個數 

	  p(0,0)=-21.0; 
	  p(0,1)=19.0; 
	  p(0,2)=-20.0;
	  p(1,0)=19.0; 
	  p(1,1)=-21.0; 
	  p(1,2)=20.0;
	  p(2,0)=40.0;
	  p(2,1)=-40.0; 
	  p(2,2)=-40.0;
}