ordinarydifferentialeguation.inl

本程序是在CCS环境下运行的

//OrdinaryDifferentialEguation.inl		求解常微分方程(组)头文件
// Ver 1.0.0.0
// 版权所有(C) 何渝, 2002
// 最后修改: 2002.5.31.

#ifndef _ORDINARYDIFFERENTIALEGUATION_INL	//避免多次编译
#define _ORDINARYDIFFERENTIALEGUATION_INL

//定步长全区间积分的欧拉法
template <class _Ty>
void ODE_EulerContentStep(_Ty t, valarray<_Ty>& y,  
										_Ty h, int k, matrix<_Ty>& z)
{
	int n = y.size();	//微分方程组中方程的个数，也是未知函数的个数
	valarray<_Ty> d(n);
    for(int i=0; i<n; i++)
		z(i,0)=y[i];
    for(int j=1; j<k; j++)
	{ 
		_Ty x=t+(j-1)*h;
        FunctionValueECS(x,y,d);
        for(i=0; i<n; i++)
          y[i]=z(i,j-1)+h*d[i];
        x=t+j*h;
        FunctionValueECS(x,y,d);
        for(i=0; i<n; i++)
          d[i]=z(i,j-1)+h*d[i];
        for(i=0; i<n; i++)
		{ 
			y[i]=(y[i]+d[i])/2.0;
            z(i,j)=y[i];
		}
	}
}

//变步长积分欧拉法
template <class _Ty>
void ODE_EulerVariationalStep(_Ty t, _Ty h, valarray<_Ty>& y, _Ty eps)
{
	int j, m(1);
    _Ty p(1.0+eps), x, q, hh=h;

	int n = y.size();	//微分方程组中方程的个数，也是未知函数的个数
    valarray<_Ty> a(n), b(n), c(n), d(n);

    for(int i=0; i<n; i++) a[i] = y[i];

    while(p>eps || FloatEqual(p,eps))
	{ 
		for(i=0; i<n; i++)
		{
			b[i]=y[i];
			y[i]=a[i];
		}
        for(j=0; j<m; j++)
		{
			for(i=0; i<n; i++) c[i]=y[i];
            x=t+j*hh;
            FunctionValueEVS(x,y,d);
            for(i=0; i<n; i++)
				y[i]=c[i]+hh*d[i];
            x=t+(j+1)*hh;
            FunctionValueEVS(x,y,d);
            for(i=0; i<n; i++)
				d[i]=c[i]+hh*d[i];
            for(i=0; i<n; i++)
				y[i]=(y[i]+d[i])/2.0;
		}
        p = 0.0;
        for(i=0; i<n; i++)
		{ 
			q=Abs(y[i]-b[i]);
            if(q>p) p=q;
        }
        hh /= 2.0;
		m = m + m;
    }
}

//全区间定步长积分维梯法
template <class _Ty>
void ODE_WittyContentStep(_Ty t, valarray<_Ty>& y,  
										_Ty h, int k, matrix<_Ty>& z)
{
    int n = y.size();	//微分方程组中方程的个数，也是未知函数的个数
	valarray<_Ty> a(n), d(n);

    for(int i=0; i<n; i++) z(i,0)=y[i];
    FunctionValueWCS(t,y,d);
    for(int j=1; j<k; j++)
    { 
		for(i=0; i<n; i++)
			a[i]=z(i,j-1)+h*d[i]/2.0;	
        _Ty x=t+(j-0.5)*h;
        FunctionValueWCS(x,a,y);
        for(i=0; i<n; i++)
        { 
			d[i]=2.0*y[i]-d[i];
            z(i,j)=z(i,j-1)+h*y[i];
        }
    }
}

//全区间定步长积分龙格-库塔法
template <class _Ty>
void ODE_RungeKuttaContentStep(_Ty t, valarray<_Ty>& y,  
										_Ty h, int k, matrix<_Ty>& z)
{
    _Ty a[4], tt;

	int n = y.size();	//微分方程组中方程的个数，也是未知函数的个数
	valarray<_Ty> b(n), d(n);

    a[0] = h / 2.0;
	a[1] = a[0];
    a[2] = a[3] = h;

    for(int i=0; i<n; i++) z(i,0)=y[i];

    for(int l=1; l<k; l++)
    { 
		FunctionValueRKCS(t, y, d);
        for(i=0; i<n; i++) b[i]=y[i];
        for(int j=0; j<=2; j++)
        {
			for(i=0; i<n; i++)
            {
				y[i]=z(i,l-1)+a[j]*d[i];
                b[i]=b[i]+a[j+1]*d[i]/3.0;
            }
            tt=t+a[j];
            FunctionValueRKCS(tt, y, d);
        }
        for(i=0; i<n; i++)
          y[i]=b[i]+h*d[i]/6.0;
        for(i=0; i<n; i++)
          z(i,l)=y[i];
        t=t+h;
    }
}

//变步长积分龙格-库塔法
template <class _Ty>
void ODE_RungeKuttaVariationalStep(_Ty t, _Ty h, 
										valarray<_Ty>& y, _Ty eps)
{
	int m(1), j, k;
    _Ty hh(h), dt,x(t),tt,q,a[4];

	int n = y.size();	//微分方程组中方程的个数，也是未知函数的个数
	valarray<_Ty> g(n), b(n), c(n), d(n), e(n);
    _Ty p=1.0+eps;

    for(int i=0; i<n; i++) c[i]=y[i];
    while(p>eps || FloatEqual(p,eps))
    {
		a[0]=hh/2.0; 
		a[1]=a[0];
		a[2] = a[3] = hh;
        for(i=0; i<n; i++)
        {
			g[i]=y[i];
			y[i]=c[i];
		}
        dt=h/m; 
		t=x;
        for(j=0; j<m; j++)
        {
			FunctionValueRKVS(t, y, d);
            for(i=0; i<n; i++) 
            {
				b[i]=y[i]; 
				e[i]=y[i];
			}
            for(k=0; k<=2; k++)
            {
				for(i=0; i<n; i++)
				{
					y[i]=e[i]+a[k]*d[i];
                    b[i]=b[i]+a[k+1]*d[i]/3.0;
                }
                tt=t+a[k];
                FunctionValueRKVS(tt, y, d);
            }
            for(i=0; i<n; i++)
              y[i]=b[i]+hh*d[i]/6.0;
            t += dt;
        }
        p = 0.0;
        for(i=0; i<n; i++)
        {
			q = Abs(y[i]-g[i]);
            if(q>p) p = q;
        }
        hh /= 2.0;
		m = m + m;
    }
}

//变步长积分基尔法
template <class _Ty>
void ODE_GillVariationalStep(_Ty t, _Ty h, valarray<_Ty>& y,
										_Ty eps, valarray<_Ty>& q)
{
    int i,k,m(1),ii;
    _Ty x(t),hh(h),r,s,t0,dt,qq;
    _Ty a[4]={0.5, 0.29289321881, 1.7071067812, 0.166666667};
    _Ty b[4]={2.0,1.0,1.0,2.0};
    _Ty c[4],e[4]={0.5,0.5,1.0,1.0};

	int n = y.size();	//微分方程组中方程的个数，也是未知函数的个数
	valarray<_Ty> g(n), u(n), v(n), d(n);

    for(i=0; i<3; i++) c[i]=a[i];
    c[3]=0.5;
    _Ty p=1.0+eps;
	
    for(int j=0; j<n; j++) u[j]=y[j];
    while(p>eps || FloatEqual(p,eps))
    {
		for(j=0; j<n; j++)
		{
			v[j]=y[j];
			y[j]=u[j];
			g[j]=q[j];
		}
        dt=h/m; 
		t=x;
        for(k=0; k<m; k++)
        { 
			FunctionValueGVS(t,y,d);
            for(ii=0; ii<=3; ii++)
            {
				for(j=0; j<n; j++)
                  d[j]=d[j]*hh;
                for(j=0; j<n; j++)
                {
					r=(a[ii]*(d[j]-b[ii]*g[j])+y[j])-y[j];
                    y[j]=y[j]+r;
                    s=g[j]+3.0*r;
                    g[j]=s-c[ii]*d[j];
                }
                t0=t+e[ii]*hh;
                FunctionValueGVS(t0,y,d);
            }
            t=t+dt;
        }
        p=0.0;
        for(j=0; j<n; j++)
        {
			qq=Abs(y[j]-v[j]);
            if(qq>p) p=qq;
        }
        hh=hh/2.0; 
		m=m+m;
    }
    for(j=0; j<n; j++) q[j]=g[j];
}

//全区间变步长积分基尔法
template <class _Ty>
void ODE_GillVariationalStepInterval(_Ty t, _Ty h, valarray<_Ty>& y,
									_Ty eps, int k, matrix<_Ty>& z)
{
    int j,m,kk,ii;
    _Ty aa(t),hh,x,p,dt,r,s,t0,qq;
    _Ty a[4]={0.5, 0.29289321881, 1.7071067812, 0.166666667};
    _Ty b[4]={2.0,1.0,1.0,2.0};
    _Ty c[4],e[4]={0.5,0.5,1.0,1.0};

	int n = y.size();	//微分方程组中方程的个数，也是未知函数的个数
	valarray<_Ty> g(n), q(n), u(n), v(n), d(n);
    
    for(int i=0; i<3; i++) c[i]=a[i];
    c[3]=0.5;
    
    for(i=0; i<n; i++) 
    {
		z(i,0) = y[i];
		q[i] = 0.0;
	}
    for(i=2; i<=k; i++)
    {
		x=aa+(i-2)*h;
		m=1;
		hh=h;
        p=1.0+eps;
        for(j=0; j<n; j++) u[j]=y[j];
        while(p>eps || FloatEqual(p,eps))
        {
			for(j=0; j<n; j++)
			{ 
				v[j]=y[j];
				y[j]=u[j];
				g[j]=q[j];
			}
            dt=h/m; 
			t=x;
            for(kk=0; kk<m; kk++)
            {
				FunctionValueGVSI(t,y,d);
                for(ii=0; ii<4; ii++)
                {
					for(j=0; j<n; j++)
                      d[j]=d[j]*hh;
                    for(j=0; j<n; j++)
                    {
						r=(a[ii]*(d[j]-b[ii]*g[j])+y[j])-y[j];
                        y[j]=y[j]+r;
                        s=g[j]+3.0*r;
                        g[j]=s-c[ii]*d[j];
                    }
                    t0=t+e[ii]*hh;
                    FunctionValueGVSI(t0,y,d);
                }
                t += dt;
            }
            p=0.0;
            for(j=0; j<n; j++)
            { 
				qq=Abs(y[j]-v[j]);
                if(qq>p) p=qq;
            }
            hh /= 2.0; 
			m=m+m;
        }
        for(j=0; j<n; j++)
        { 
			q[j]=g[j];
			z(j,i-1)=y[j];
		}
    }
}

//变步长积分默森法
template <class _Ty>
void ODE_MersonVariationalStep(_Ty t, _Ty h, int n, valarray<_Ty>& y,
										_Ty eps, int k, matrix<_Ty>& z)
{
    int j,m,nn;
    _Ty aa(t),bb,x,hh,p,dt,t0,qq;
    valarray<_Ty> a(n), b(n), c(n), u(n), v(n), d(n);
        
    for(int i=0; i<n; i++) z(i,0)=y[i];
    for(i=1; i<k; i++)
    { 
		x=aa+(i-1)*h; 
		nn=1; 
		hh=h;
        for(j=0; j<n; j++) u[j]=y[j];
        p=1.0+eps;
        while(p>eps || FloatEqual(p,eps))
        {
			for(j=0; j<n; j++)
			{
				v[j]=y[j]; 
				y[j]=u[j];
			}
            dt=h/nn; t=x;
            for(m=0; m<nn; m++)
			{ 
				FunctionValueMVS(t,y,d);
                for(j=0; j<n; j++)
                {
					a[j]=d[j];
					y[j]=y[j]+hh*d[j]/3.0;
				}
                t0=t+hh/3.0;
                FunctionValueMVS(t0,y,d);
                for(j=0; j<n; j++)
                {
					b[j]=d[j];
					y[j]=y[j]+hh*(d[j]-a[j])/6.0;
				}
                FunctionValueMVS(t0,y,d);
                for(j=0; j<n; j++)
                {
					b[j]=d[j];
                    bb=(d[j]-4.0*(b[j]+a[j]/4.0)/9.0)/8.0;
                    y[j]=y[j]+3.0*hh*bb;
                }
                t0=t+hh/2.0;
                FunctionValueMVS(t0,y,d);
                for(j=0; j<n; j++)
                { c[j]=d[j];
                    qq=d[j]-15.0*(b[j]-a[j]/5.0)/16.0;
                    y[j]=y[j]+2.0*hh*qq;
                }
                t0=t+hh;
                FunctionValueMVS(t0,y,d);
                for(j=0; j<n; j++)
                { 
					qq=c[j]-9.0*(b[j]-2.0*a[j]/9.0)/8.0;
                    qq=d[j]-8.0*qq;
                    y[j]=y[j]+hh*qq/6.0;
                }
                t += dt;
            }
            p=0.0;
            for(j=0; j<n; j++)
            { qq=Abs(y[j]-v[j]);
                if(qq>p) p=qq;
            }
            hh /= 2.0;
			nn=nn+nn;
        }
        for(j=0; j<n; j++) z(j,i)=y[j];
    }
}

//积分一步连分式法
template <class _Ty>
void ODE_Fraction(_Ty t, _Ty h, valarray<_Ty>& y, _Ty eps)
{
	int i,j,k(1),m,nn(1),it(1);

	int n = y.size();	//微分方程组中方程的个数，也是未知函数的个数
    _Ty x(t),hh(h),dd,q,p,g[10];
	valarray<_Ty> e(n), b(10*n), w(n), u(n), v(n), d(n);

    for(j=0; j<n; j++) v[j]=y[j];
	g[0]=hh;
    Fraction1(x,hh,y,w,d,e);
    for(j=0; j<n; j++)
    { 
		b[j]=y[j];
		u[j]=y[j];
	}

    while(it==1)
    { 
		nn=nn+nn; 
		hh=hh/2.0;
		it=0;
        g[k]=hh;
        for(j=0; j<n; j++) y[j]=v[j];
        t=x;
        for(j=0; j<nn; j++)
        { 
			Fraction1(t,hh,y,w,d,e);
            t=t+hh;
        }
        for(j=0; j<n; j++)
        { 
			dd=y[j]; 
			m=0;
            for(i=0; i<k; i++)
              if(m==0)
              {
				  q=dd-b[i*n+j];
                  if(FloatEqual(q,0.0)) m=1;
                  else dd=(g[k]-g[i])/q;
              }
            b[k*n+j]=dd;
            if(m!=0) b[k*n+j]=1.0e+35;
        }
        for(j=0; j<n; j++)
        {
			dd=0.0;
            for(i=k-1; i>=0; i--)
				dd=-g[i]/(b[(i+1)*n+j]+dd);
            y[j]=dd+b[j];
        }
        p=0.0;
        for(j=0; j<n; j++)
        { 
			q=Abs(y[j]-u[j]);
            if(q>p) p=q;
        }
        if((p>eps || FloatEqual(p,eps))&&(k<7))
        { 
			for(j=0; j<n; j++) u[j]=y[j];
            k=k+1; 
			it=1;
        }
    }
}

//积分一步连分式法ODE_Fraction()辅助函数
template <class _Ty>
void Fraction1(_Ty t, _Ty h, valarray<_Ty>& y, valarray<_Ty>& b, 
									valarray<_Ty>& d, valarray<_Ty>& e)
{
    _Ty a[4],tt;
	int n = y.size();	//微分方程组中方程的个数，也是未知函数的个数
    
	a[0]=h/2.0;
	a[1]=a[0];
	a[2]=h; 
	a[3]=h;
    FunctionValueF(t,y,d);
    for(int i=0; i<n; i++)
	{
		b[i]=y[i];
		e[i]=y[i];
	}
    for(int k=0; k<=2; k++)
    {
		for(i=0; i<n; i++)
        {
			y[i]=e[i]+a[k]*d[i];
            b[i]=b[i]+a[k+1]*d[i]/3.0;
        }
        tt=t+a[k];
        FunctionValueF(tt,y,d);
    }
    for(i=0; i<n; i++)
      y[i]=b[i]+h*d[i]/6.0;
}


//全区间积分连分式法
template <class _Ty>
void ODE_FractionInterval(_Ty t, _Ty h, valarray<_Ty>& y,
									_Ty eps, int k, matrix<_Ty>& z)
{
	int j,kk,l,m,nn,it;
    _Ty x,hh,dd,tt,q,p,g[10];

	int n = y.size();	//微分方程组中方程的个数，也是未知函数的个数
	valarray<_Ty> e(n), b(10*n), w(n), u(n), v(n), d(n);

    for(int i=0; i<n; i++) z(i,0)=y[i];
    for(i=1; i<k; i++)
    { 
		for(j=0; j<n; j++) v[j]=y[j];
        x=t+(i-1)*h; 
		nn=1; 
		hh=h; 
		g[0]=hh;
        Fraction2(x,hh,y,w,d,e);
        for(j=0; j<n; j++)