ordinarydifferentialeguation.inl

矩阵运算

        {
			b[j]=y[j];
			u[j]=y[j];
		}
        kk=1;
		it=1;
        while(it==1)
        {
			nn=nn+nn; 
			hh=hh/2.0;
			it=0;
            g[kk]=hh;
            for(j=0; j<n; j++) y[j]=v[j];
            tt=x;
            for(j=0; j<nn; j++)
            {
				Fraction2(t,hh,y,w,d,e);
                tt=tt+hh;
            }
            for(j=0; j<n; j++)
            { 
				dd=y[j]; 
				l=0;
                for(m=0; m<kk; m++)
                  if(l==0)
                  {
					  q=dd-b[m*n+j];
                      if(Abs(q)+1.0==1.0) l=1;
                      else dd=(g[kk]-g[m])/q;
                  }
                b[kk*n+j]=dd;
                if(l!=0) b[kk*n+j]=1.0e+35;
            }
            for(j=0; j<n; j++)
            { 
				dd=0.0;
                for(m=kk-1; m>=0; m--)
					dd=-g[m]/(b[(m+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))&&(kk<7))
            { for(j=0; j<n; j++) u[j]=y[j];
                kk=kk+1; it=1;
            }
        }
        for(j=0; j<n; j++)
          z(j,i)=y[j];
    }
}

//全区间积分连分式法辅助函数
template <class _Ty>
void Fraction2(_Ty t, _Ty h, valarray<_Ty>& y, valarray<_Ty>& b, 
									valarray<_Ty>& d, valarray<_Ty>& e)
{
    int i,k;
    _Ty a[4],tt;

	int n = y.size();	//微分方程组中方程的个数，也是未知函数的个数
    a[0]=h/2.0; 
	a[1]=a[0];
	a[2]=h;
	a[3]=h;
    FunctionValueFI(t,y,d);
    for(i=0; i<n; i++) 
	{
		b[i]=y[i]; 
		e[i]=y[i];
	}
    for(k=0; k<3; 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];
        FunctionValueFI(tt,y,d);
    }
    for(i=0; i<n; i++)
      y[i]=b[i]+h*d[i]/6.0;
}

//全区间积分双边法
template <class _Ty>
void ODE_BothSidesInterval(_Ty t, _Ty h, valarray<_Ty>& y,
					   				_Ty eps, int k, matrix<_Ty>& z)
{
    _Ty a(t), qq;

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

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

    FunctionValueBSI(t,y,d);
    for(int j=0; j<n; j++) u[j]=d[j];
    BothSidesInterval(t,h,y,eps);		//辅助函数
    t=a+h;
    FunctionValueBSI(t,y,d);
    for(j=0; j<n; j++)
    { 
		z(j,1)=y[j]; v[j]=d[j];}
    for(j=0; j<n; j++)
    { 
		p[j]=-4.0*z(j,1)+5.0*z(j,0)+2.0*h*(2.0*v[j]+u[j]);
        y[j]=p[j];
    }
    t=a+2.0*h;
    FunctionValueBSI(t,y,d);
    for(j=0; j<n; j++)
    { 
		qq=2.0*h*(d[j]-2.0*v[j]-2.0*u[j])/3.0;
        qq=qq+4.0*z(j,1)-3.0*z(j,0);
        z(j,2)=(p[j]+qq)/2.0;
        y[j]=z(j,2);
    }
    for(i=3; i<k; i++)
    { 
		t=a+(i-1)*h;
        FunctionValueBSI(t,y,d);
        for(j=0; j<n; j++)
        {
			u[j]=v[j];
			v[j]=d[j];
		}
        for(j=0; j<n; j++)
        { 
			qq=-4.0*z(j,i-1)+5.0*z(j,i-2);
            p[j]=qq+2.0*h*(2.0*v[j]+u[j]);
            y[j]=p[j];
        }
        t=t+h;
        FunctionValueBSI(t,y,d);
        for(j=0; j<n; j++)
        { 
			qq=2.0*h*(d[j]-2.0*v[j]-2.0*u[j])/3.0;
            qq=qq+4.0*z(j,i-1)-3.0*z(j,i-2);
            y[j]=(p[j]+qq)/2.0;
            z(j,i)=y[j];
        }
    }
}

//全区间积分双边法辅助函数
template <class _Ty>
void BothSidesInterval(_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> d(n), b(n), c(n), e(n), g(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]=hh;
		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++)
        {
			FunctionValueBSI(t,y,d);
            for(i=0; i<n; i++) 
            {
				b[i]=y[i];
				e[i]=y[i];
			}
            for(k=0; k<3; 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];
                FunctionValueBSI(tt,y,d);
            }
            for(i=0; i<n; i++)
              y[i]=b[i]+hh*d[i]/6.0;
            t=t+dt;
        }
        p=0.0;
        for(i=0; i<n; i++)
        { 
			q=Abs(y[i]-g[i]);
            if(q>p) p=q;
        }
        hh=hh/2.0;
		m=m+m;
    }
}

//全区间积分阿当姆斯预报-校正法
template <class _Ty>
void ODE_AdamsInterval(_Ty t, _Ty h, valarray<_Ty>& y,
									_Ty eps, int k, matrix<_Ty>& z)
{
    int j,m;
    _Ty a(t), q;

	int n = y.size();	//微分方程组中方程的个数，也是未知函数的个数
	valarray<_Ty> b(4*n), e(n), s(n), g(n), d(n);
    
    for(int i=0; i<n; i++) z(i,0)=y[i];
    FunctionValueAI(t,y,d);
    for(i=0; i<n; i++) b[i]=d[i];
    for(i=1; i<4; i++)
      if(i<k)
      { 
		  t=a+i*h;
          AdamsInterval(t,h,y,eps);			//辅助函数
          for(j=0; j<n; j++) z(j,i)=y[j];
          FunctionValueAI(t,y,d);
          for(j=0; j<n; j++) b[i*n+j]=d[j];
      }
    for(i=4; i<k; i++)
    { 
		for(j=0; j<n; j++)
        { 
			q=55.0*b[3*n+j]-59.0*b[2*n+j];
            q=q+37.0*b[n+j]-9.0*b[j];
            y[j]=z(j,i-1)+h*q/24.0;
            b[j]=b[n+j];
            b[n+j]=b[n+n+j];
            b[n+n+j]=b[n+n+n+j];
        }
        t=a+i*h;
        FunctionValueAI(t,y,d);
        for(m=0; m<n; m++) b[n+n+n+m]=d[m];
        for(j=0; j<n; j++)
        { 
			q=9.0*b[3*n+j]+19.0*b[n+n+j]-5.0*b[n+j]+b[j];
            y[j]=z(j,i-1)+h*q/24.0;
            z(j,i)=y[j];
        }
        FunctionValueAI(t,y,d);
        for(m=0; m<n; m++) b[3*n+m]=d[m];
    }
}

//全区间积分阿当姆斯预报-校正法辅助函数
template <class _Ty>
void AdamsInterval(_Ty t, _Ty h, valarray<_Ty>& y, _Ty eps)
{
    int m(1), j, k;

	int n = y.size();	//微分方程组中方程的个数，也是未知函数的个数
    _Ty hh(h), dt, x(t), tt, q, a[4];
	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]=hh;
		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++)
        { 
			FunctionValueAI(t,y,d);
            for(i=0; i<n; i++) 
            {
				b[i]=y[i];
				e[i]=y[i];
			}
            for(k=0; k<3; 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];
                FunctionValueAI(tt,y,d);
            }
            for(i=0; i<n; i++)
              y[i]=b[i]+hh*d[i]/6.0;
            t=t+dt;
        }
        p=0.0;
        for(i=0; i<n; i++)
        {
			q=Abs(y[i]-g[i]);
            if(q>p) p=q;
        }
        hh=hh/2.0;
		m=m+m;
    }
}

//全区间积分哈明法
template <class _Ty>
void ODE_HammingInterval(_Ty t, _Ty h, valarray<_Ty>& y,
										_Ty eps, int k, matrix<_Ty>& z)
{
    int j,m;
    _Ty a(t), q;

	int n = y.size();	//微分方程组中方程的个数，也是未知函数的个数
	valarray<_Ty> d(n), w(n), u(n), v(n), b(4*n), g(n);
    
    for(int i=0; i<n; i++) z(i,0)=y[i];
    FunctionValueHI(t,y,d);
    for(i=0; i<n; i++) b[i]=d[i];
    for(i=1; i<=3; i++)
      if(i<k)
      { 
		  t=a+i*h;
          HammingInterval(t,h,y,eps);
          for(m=0; m<n; m++) z(m,i)=y[m];
          FunctionValueHI(t,y,d);
          for(m=0; m<n; m++) b[i*n+m]=d[m];
      }
    for(i=0; i<n; i++) u[i]=0.0;
    for(i=4; i<k; i++)
    { 
		for(j=0; j<n; j++)
        {
			q=2.0*b[3*n+j]-b[n+n+j]+2.0*b[n+j];
            y[j]=z(j,i-4)+4.0*h*q/3.0;
        }
        for(j=0; j<n; j++)
			y[j]=y[j]+112.0*u[j]/121.0;
        t=a+i*h;
        FunctionValueHI(t,y,d);
        for(j=0; j<n; j++)
        { 
			q=9.0*z(j,i-1)-z(j,i-3);
            q=(q+3.0*h*(d[j]+2.0*b[3*n+j]-b[n+n+j]))/8.0;
            u[j]=q-y[j];
            z(j,i)=q-9.0*u[j]/121.0;
            y[j]=z(j,i);
            b[n+j]=b[n+n+j];
            b[n+n+j]=b[n+n+n+j];
        }
        FunctionValueHI(t,y,d);
        for(m=0; m<n; m++) b[3*n+m]=d[m];
    }
}

//全区间积分哈明法辅助函数
template <class _Ty>
void HammingInterval(_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]=hh;
		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++)
        {
			FunctionValueHI(t,y,d);
            for(i=0; i<n; i++) 
            {
				b[i]=y[i];
				e[i]=y[i];
			}
            for(k=0; k<3; 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];
                FunctionValueHI(t,y,d);
            }
            for(i=0; i<n; i++)
				y[i]=b[i]+hh*d[i]/6.0;
            t=t+dt;
        }
        p=0.0;
        for(i=0; i<n; i++)
        { 
			q=Abs(y[i]-g[i]);
            if(q>p) p=q;
        }
        hh=hh/2.0; 
		m=m+m;
    }
}

//积分一步特雷纳法
template <class _Ty>
void ODE_Treanor(_Ty t, _Ty h, valarray<_Ty>& y)
{
//extern void gtnr1f();
    _Ty s,aa,bb,dd,g,dy,dy1;

	int n = y.size();	//微分方程组中方程的个数，也是未知函数的个数
	valarray<_Ty> w(4*n), q(4*n), r(4*n), d(n), p(n);
    
    for(int j=0; j<n; j++) w[j]=y[j];
    FunctionValueT(t,y,d);
    for(j=0; j<n; j++)
    {
		q[j]=d[j]; 
		y[j]=w[j]+h*d[j]/2.0;
        w[n+j]=y[j];
    }
    s=t+h/2.0;
    FunctionValueT(t,y,d);
    for(j=0; j<n; j++)
    { 
		q[n+j]=d[j];
        y[j]=w[j]+h*d[j]/2.0;
        w[n+n+j]=y[j];
    }
    FunctionValueT(t,y,d);
    for(j=0; j<n; j++) q[n+n+j]=d[j];

    for(j=0; j<n; j++)
    { 
		aa=q[n+n+j]-q[n+j];
        bb=w[n+n+j]-w[n+j];
        if(-aa*bb*h>0.0)
        {
			p[j]=-aa/bb; dd=-p[j]*h;
            r[j]=exp(dd);
            r[n+j]=(r[j]-1.0)/dd;
            r[n+n+j]=(r[n+j]-1.0)/dd;
            r[3*n+j]=(r[n+n+j]-1.0)/dd;
        }
        else p[j]=0.0;
        if(p[j]<0.0||FloatEqual(p[j],0.0)) g=q[n+n+j];
        else
        { 
			g=2.0*(q[n+n+j]-q[j])*r[n+n+j];
            g=g+(q[j]-q[n+j])*r[n+j]+q[n+j];
        }
        w[3*n+j]=w[j]+g*h;
        y[j]=w[3*n+j];
    }
    s=t+h;
    FunctionValueT(t,y,d);
    for(j=0; j<n; j++) q[3*n+j]=d[j];

    for(j=0; j<n; j++)
    { 
		if(p[j]<0.0||FloatEqual(p[j],0.0))
        { 
			dy=q[j]+2.0*(q[n+j]+q[n+n+j]);
            dy=(dy+q[n+n+n+j])*h/6.0;
        }
        else
        {
			dy=-3.0*(q[j]+p[j]*w[j])+2.0*(q[n+j]+p[j]*w[n+j]);
            dy=dy+2.0*(q[n+n+j]+p[j]*w[n+n+j]);
            dy=dy-(q[n+n+n+j]+p[j]*w[n+n+n+j]);
            dy=dy*r[n+n+j]+q[j]*r[n+j];
            dy1=q[j]-q[n+j]-q[n+n+j]+q[n+n+n+j];
            dy1=dy1+(w[j]-w[n+j]-w[n+n+j]+w[n+n+n+j])*p[j];
            dy=(dy+4.0*dy1*r[n+n+n+j])*h;
        }
        y[j]=w[j]+dy;
    }
}

//全区间积分特雷纳法
template <class _Ty>
void ODE_TreanorInterval(_Ty t, _Ty h, valarray<_Ty>& y,
											int k, matrix<_Ty>& z)
{
    _Ty a(t),s,aa,bb,dd,g,dy,dy1;

	int n = y.size();	//微分方程组中方程的个数，也是未知函数的个数
	valarray<_Ty> w(4*n), q(4*n), r(4*n), d(n), p(n);
    
    for(int j=0; j<n; j++) z(j,0)=y[j];

    for(int i=1; i<k; i++)
    {
		t=a+(i-1)*h;
        for(j=0; j<n; j++) w[j]=y[j];
        FunctionValueTI(s,y,d);
        for(j=0; j<n; j++)
        { 
			q[j]=d[j]; 
			y[j]=w[j]+h*d[j]/2.0;
            w[n+j]=y[j];
        }
        s=t+h/2.0;
        FunctionValueTI(s,y,d);
        for(j=0; j<n; j++)
        { 
			q[n+j]=d[j];
            y[j]=w[j]+h*d[j]/2.0;
            w[n+n+j]=y[j];
        }