rk6.cpp

来自「使用runge-kutta方法求解ODE问题」· C++ 代码 · 共 316 行

#include <iostream>
#include <fstream>
#include <math.h>
#include <cstdlib>
#include <vector>
#define N 8000
#define PI 4*atan(1.0)
#define h (2*PI)/N

using namespace std;

 double y[N+1][2];   //这是用来存储N+1层的数值解
double t[N+1];          //这是用来存储N+1层的时间
 double y_exact[N+1][2]; //这是用来存储N+1层的真解
void gausssseidel0( double aa[6][6], double x[6],double bb[6]); // 这是求解二阶的线性方程组的函数
void gausssseidel( double aa[4][4], double x[4], double bb[4]); // 这是求解四阶的线性代数方程组的函数
//下面是算出矩阵J的四个元素，变量是时间 t
double J00(double t); 
double J01(double t);
double J10(double t);
double J11(double t);
//下面是算出向量g的四个元素，变量是时间 t
double g0(double t);
double g1(double t);


int main()
{
	int i;
	 double cof[6][6];
	 double cof1[4][4];
           double rhs1[4];
	  double rhs[6];
         double Y[4];
        double c00=0.0, c01=0.0, c10=5.0/4, c11=1.0/4; //这是矩阵C11的四个元素
        double d00=0.0, d01=1.0, d10=1.0/2, d11=1.0/2;  //这是矩阵C12的四个元素
        double e00=0.0, e01=0.0, e10=4.0/5, e11=2.0/5; //这是矩阵C21的四个元素
        double f00=0.0, f01=1.0, f10=1.0/5, f11=4.0/5; //这是矩阵C22的四个元素
     //下面是初始解
	y[0][0]=1;
	y[0][1]=0;
     //下面是0时刻的真解
	y_exact[0][0]=1;
	y_exact[0][1]=0;
	t[0]=0.0;
    //下面是用三级Runge-Kutta方法求第一时间层的数值解的那些系数
       double a11=1.0/6;
       double a12=-1.0/3;
       double a13=1.0/6;
       double a21=1.0/6;
       double a22=5.0/12;
       double a23=-1.0/12;
       double a31=1.0/6;
       double a32=2.0/3;
       double a33=1.0/6;
       double c1=0;
       double c2=1.0/2;
       double c3=1.0;
       double b2=2.0/3;
       double b3=1.0/6;
       double b1=1.0/6;
       //下面就是按照三级Runge-Kutta法求出第一时间层的数值解。
       double  Y1[2], Y2[2], Y3[2],Y4[6];

       cof[0][0]=1-h*a11*J00(t[0]+c1*h);
       cof[0][1]=-h*a11*J01(t[0]+c1*h);
       cof[0][2]=-h*a12*J00(t[0]+c2*h);
       cof[0][3]=-h*a12*J01(t[0]+c2*h);
       cof[0][4]=-h*a13*J00(t[0]+c3*h);
       cof[0][5]=-h*a13*J01(t[0]+c3*h);

       cof[1][0]=-h*a11*J10(t[0]+c1*h);
       cof[1][1]=1-h*a11*J11(t[0]+c1*h);
       cof[1][2]=-h*a12*J10(t[0]+c2*h);
       cof[1][3]=-h*a12*J11(t[0]+c2*h);
       cof[1][4]=-h*a13*J10(t[0]+c3*h);
       cof[1][5]=-h*a13*J11(t[0]+c3*h);

       cof[2][0]=-h*a21*J00(t[0]+c1*h);
       cof[2][1]=-h*a21*J01(t[0]+c1*h);
       cof[2][2]=1-h*a22*J00(t[0]+c2*h);
       cof[2][3]=-h*a22*J01(t[0]+c2*h);
       cof[2][4]=-h*a23*J00(t[0]+c3*h);
       cof[2][5]=-h*a23*J01(t[0]+c3*h);

       cof[3][0]=-h*a21*J10(t[0]+c1*h);
       cof[3][1]=-h*a21*J11(t[0]+c1*h);
       cof[3][2]=-h*a22*J10(t[0]+c2*h);
       cof[3][3]=1-h*a22*J11(t[0]+c2*h);
       cof[3][4]=-h*a23*J10(t[0]+c3*h);
       cof[3][5]=-h*a23*J11(t[0]+c3*h);

       cof[4][0]=-h*a31*J00(t[0]+c1*h);
       cof[4][1]=-h*a31*J01(t[0]+c1*h);
       cof[4][2]=-h*a32*J00(t[0]+c2*h);
       cof[4][3]=-h*a32*J01(t[0]+c2*h);
       cof[4][4]=1-h*a33*J00(t[0]+c3*h);
       cof[4][5]=-h*a33*J01(t[0]+c3*h);

       cof[5][0]=-h*a31*J10(t[0]+c1*h);
       cof[5][1]=-h*a31*J11(t[0]+c1*h);
       cof[5][2]=-h*a32*J10(t[0]+c2*h);
       cof[5][3]=-h*a32*J11(t[0]+c2*h);
       cof[5][4]=-h*a33*J10(t[0]+c3*h);
       cof[5][5]=1-h*a33*J11(t[0]+c3*h);

       
       rhs[0]=y[0][0]+h*a11*g0(t[0]+c1*h)+h*a12*g0(t[0]+c2*h)+h*a13*g0(t[0]+c3*h);
       rhs[1]=y[0][1]+h*a11*g1(t[0]+c1*h)+h*a12*g1(t[0]+c2*h)+h*a13*g1(t[0]+c3*h);;
       
       rhs[2]=y[0][0]+h*a21*g0(t[0]+c1*h)+h*a22*g0(t[0]+c2*h)+h*a23*g0(t[0]+c3*h);
       rhs[3]=y[0][1]+h*a21*g1(t[0]+c1*h)+h*a22*g1(t[0]+c2*h)+h*a23*g1(t[0]+c3*h);
       
       rhs[4]=y[0][0]+h*a31*g0(t[0]+c1*h)+h*a32*g0(t[0]+c2*h)+h*a33*g0(t[0]+c3*h);
       rhs[5]=y[0][1]+h*a31*g1(t[0]+c1*h)+h*a32*g1(t[0]+c2*h)+h*a33*g1(t[0]+c3*h);

       gausssseidel0(cof,Y4, rhs);

       Y1[0]=Y4[0];
       Y1[1]=Y4[1];
       Y2[0]=Y4[2];
       Y2[1]=Y4[3];
       Y3[0]=Y4[4];
       Y3[1]=Y4[5];

       y[1][0]=y[0][0]+h*b1*J00(t[0]+c1*h)*Y1[0]+h*b1*J01(t[0]+c1*h)*Y1[1]+h*b2*J00(t[0]+c2*h)*Y2[0]+h*b2*J01(t[0]+c2*h)*Y2[1]+h*b3*J00(t[0]+c3*h)*Y3[0]+h*b3*J01(t[0]+c3*h)*Y3[1]+h*b1*g0(t[0]+c1*h)+h*b2*g0(t[0]+c2*h)+h*b3*g0(t[0]+c3*h);
       y[1][1]=y[0][1]+h*b1*J10(t[0]+c1*h)*Y1[0]+h*b1*J11(t[0]+c1*h)*Y1[1]+h*b2*J10(t[0]+c2*h)*Y2[0]+h*b2*J11(t[0]+c2*h)*Y2[1]+h*b3*J10(t[0]+c3*h)*Y3[0]+h*b3*J11(t[0]+c3*h)*Y3[1]+h*b1*g1(t[0]+c1*h)+h*b2*g1(t[0]+c2*h)+h*b3*g1(t[0]+c3*h);   
           //计算第一时间层的真解
        t[1]=t[0]+h;  
	y_exact[1][0]=cos(t[1]);
	y_exact[1][1]=sin(t[1]);
/*
      //计算第一时间层真解与数值解的绝对误差
	double error00=0.0,error10=0.0;
	for(i=1;i<=1;i++)
	{
		double e00=0.0,e10=0.0;
		e00=(y_exact[i][0]-y[i][0])*(y_exact[i][0]-y[i][0]);
		e10=(y_exact[i][1]-y[i][1])*(y_exact[i][1]-y[i][1]);
		error00=error00+e00;
		error10=error10+e10;
	}
	cout<<"error00="<<sqrt(error00+error10)/2<<endl;

    */     
        double mu1=1.0;
        double mu2=2.0;
            //下面就是根据2级2步方法给出的方程求出后面时间层的数值解。
	for (i=2; i<=N; i++)
	{
		t[i]=t[i-1]+h;
               //计算利用公式1求解Y的线性代数方程组的系数。
		cof1[0][0]=1;
		cof1[0][1]=0;
		cof1[0][2]=0;
		cof1[0][3]=0;
		cof1[1][0]=0;
		cof1[1][1]=1;
		cof1[1][2]=0;
		cof1[1][3]=0;
		cof1[2][0]=-h*c10*J00(t[i-2]+mu1*h);
		cof1[2][1]=-h*c10*J01(t[i-2]+mu1*h);
		cof1[2][2]=1-h*(c11*J00(t[i-2]+mu2*h));
		cof1[2][3]=-h*(c11*J01(t[i-2]+mu2*h));
		cof1[3][0]=-h*c10*J10(t[i-2]+mu1*h);
		cof1[3][1]=-h*c10*J11(t[i-2]+mu1*h);
		cof1[3][2]=-h*c11*J10(t[i-2]+mu2*h);
		cof1[3][3]=1-h*c11*J11(t[i-2]+mu2*h);
            //计算右端项
                rhs1[0]=h*(c00*g0(t[i-2]+mu1*h)+c01*g0(t[i-2]+mu2*h))+d00*y[i-2][0]+d01*y[i-1][0];
                rhs1[1]=h*(c00*g1(t[i-2]+mu1*h)+c01*g1(t[i-2]+mu2*h))+d00*y[i-2][1]+d01*y[i-1][1];
                rhs1[2]=h*(c10*g0(t[i-2]+mu1*h)+c11*g0(t[i-2]+mu2*h))+d10*y[i-2][0]+d11*y[i-1][0];
                rhs1[3]=h*(c10*g1(t[i-2]+mu1*h)+c11*g1(t[i-2]+mu2*h))+d10*y[i-2][1]+d11*y[i-1][1];
           // 调用高斯塞德尔迭代法求解方程
		gausssseidel(cof1,Y, rhs1);

               //利用上面求得的Y，根据第二个公式求出y.
                y[i][0]=f10*y[i-2][0]+f11*y[i-1][0]+h*(e10*J00(t[i-2]+mu1*h)*Y[0]+e10*J01(t[i-2]+mu1*h)*Y[1]+e11*J00(t[i-2]+mu2*h)*Y[2]+e11*J01(t[i-2]+mu2*h)*Y[3])+h*(e10*g0(t[i-2]+mu1*h)+e11*g0(t[i-2]+mu2*h));
                y[i][1]=f10*y[i-2][1]+f11*y[i-1][1]+h*(e10*J10(t[i-2]+mu1*h)*Y[0]+e10*J11(t[i-2]+mu1*h)*Y[1]+e11*J10(t[i-2]+mu2*h)*Y[2]+e11*J11(t[i-2]+mu2*h)*Y[3])+h*(e10*g1(t[i-2]+mu1*h)+e11*g1(t[i-2]+mu2*h));

		y_exact[i][0]=cos(t[i]);
		y_exact[i][1]=sin(t[i]);
	}
	// 计算所有时间层的误差
	double error0=0.0,error1=0.0;
        ofstream os;
	os.open("l2error8000.txt");
        os<<"The discrete l2-norm error on each step are :"<<endl;
         os<<"The time step h :"<<"\t\t"<<h<<"\t\t"<<"The N :"<<"\t\t"<<N<<"\n";
	for(i=0;i<=N;i++)
	{
		double e0,e1;
		e0=(y_exact[i][0]-y[i][0])*(y_exact[i][0]-y[i][0]);
		e1=(y_exact[i][1]-y[i][1])*(y_exact[i][1]-y[i][1]);
               os<<"t["<<i<<"]"<<"\t\t"<<sqrt(e0+e1)<<"\n";
		error0=error0+e0;
		error1=error1+e1;
	}
         os<<"The total l2 error is :"<<"\t\t"<<sqrt(error0+error1)<<"\n";
       os.close();
      cout<<"h="<<h<<endl;
	cout<<"error0="<<sqrt(error0)<<"\t"<<"error1="<<sqrt(error1)<<endl;
/*
        double error[4];
	for(i=0;i<=4;i++)
	{
               error[i]=0.0;
         }
		error[0]+=(y_exact[1250][0]-y[1250][0])*(y_exact[1250][0]-y[1250][0])+(y_exact[1250][1]-y[1250][1])*(y_exact[1250][1]-y[1250][1]);
		error[1]+=(y_exact[2500][0]-y[2500][0])*(y_exact[2500][0]-y[2500][0])+(y_exact[2500][1]-y[2500][1])*(y_exact[2500][1]-y[2500][1]);
                error[2]+=(y_exact[3750][0]-y[3750][0])*(y_exact[3750][0]-y[3750][0])+(y_exact[3750][1]-y[3750][1])*(y_exact[3750][1]-y[3750][1]);
                error[3]+=(y_exact[5000][0]-y[5000][0])*(y_exact[5000][0]-y[5000][0])+(y_exact[5000][1]-y[5000][1])*(y_exact[5000][1]-y[5000][1]);
	cout<<"error2="<<sqrt(error[0])/2<<"\t"<<"error3="<<sqrt(error[1])/2<<endl;
       cout<<"error4="<<sqrt(error[2])/2<<"\t"<<"error5="<<sqrt(error[3])/2<<endl;

               double error5=0.0;
              double error6=0.0;
                error5+=(y_exact[250][0]-y[250][0])*(y_exact[250][0]-y[250][0])+(y_exact[250][1]-y[250][1])*(y_exact[250][1]-y[250][1]);
                error6+=(y_exact[500][0]-y[500][0])*(y_exact[500][0]-y[500][0])+(y_exact[500][1]-y[500][1])*(y_exact[500][1]-y[500][1]);
	cout<<"error5="<<sqrt(error5)/2<<"\t"<<"error6="<<sqrt(error6)/2<<endl;
*/

   return 1;
}

//高斯塞德尔迭代方法求二阶方程组。
void gausssseidel0(double aa[6][6],double x[6],double bb[6])
{
    double error;
   double  b,c,d;
   int i,j;  
   do{
    error=0.0; 
   for (i=0; i<6; i++)
   {
        c=x[i];
      for(j=0;j<6;j++)
      {
        if (j!=i)
        {
              b+=aa[i][j]*x[j];
        }   
     }    
         x[i]=(bb[i]-b)/aa[i][i];
         b=0;
          d=x[i];
       if(fabs(c-d)>error)
                           error=fabs(c-d);
   }
 //        cout<<"err="<<error<<endl;
    if (error<1.0e-08) break;
   }while(1);
}


//高斯塞德尔迭代方法求四阶方程组。
void gausssseidel(double aa[4][4], double x[4],double bb[4])
{
    double error;
   double  b,c,d;
   int i,j;  
   do{
    error=0.0; 
   for (i=0; i<4; i++)
   {
        c=x[i];
      for(j=0;j<4;j++)
      {
        if (j!=i)
        {
              b+=aa[i][j]*x[j];
        }   
     }    
         x[i]=(bb[i]-b)/aa[i][i];
         b=0;
          d=x[i];
       if(fabs(c-d)>error)
                           error=fabs(c-d);
   }
 //        cout<<"err="<<error<<endl;
    if (error<1.0e-08) break;
   }while(1);
}


double J00(double t)
{
       return -1.0+sin(exp(-t));
}

double J01(double t)
{
       return 0.0;
}

double J10(double t)
{
       return 10;
}

double J11(double t)
{
       return -1.0e3*exp(-t/100)+sin(exp(-t));
}

double g0(double t)
{
       return cos(t)*(1-sin(exp(-t)))-sin(t);
}

double g1(double t)
{
       return  (1000*exp(-t/100)-sin(exp(-t)))*sin(t)-9*cos(t);
}