开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 数值算法/人工智能 > 查看源码
OrdinaryDifferentialEguation.inl
资源名称：C++常用算法之06矩阵特征值与特征向量的计算.rar　[点击查看]
上传用户：fxromeo
上传日期：2010-04-08
资源大小：89k
文件大小：36k
源码类别：

数值算法/人工智能

开发平台：

Visual C++

OrdinaryDifferentialEguation.inl：源码内容
  1. //OrdinaryDifferentialEguation.inl 求解常微分方程(组)头文件
  2. // Ver 1.0.0.0
  3. // 版权所有(C) 何渝, 2002
  4. // 最后修改: 2002.5.31.
  6. #ifndef _ORDINARYDIFFERENTIALEGUATION_INL //避免多次编译
  7. #define _ORDINARYDIFFERENTIALEGUATION_INL
  9. //定步长全区间积分的欧拉法
  10. template <class _Ty>
  11. void ODE_EulerContentStep(_Ty t, valarray<_Ty>& y,  
  12. _Ty h, int k, matrix<_Ty>& z)
  13. {
  14. int n = y.size(); //微分方程组中方程的个数，也是未知函数的个数
  15. valarray<_Ty> d(n);
  16.     for(int i=0; i<n; i++)
  17. z(i,0)=y[i];
  18.     for(int j=1; j<k; j++)
  19. { 
  20. _Ty x=t+(j-1)*h;
  21.         FunctionValueECS(x,y,d);
  22.         for(i=0; i<n; i++)
  23.           y[i]=z(i,j-1)+h*d[i];
  24.         x=t+j*h;
  25.         FunctionValueECS(x,y,d);
  26.         for(i=0; i<n; i++)
  27.           d[i]=z(i,j-1)+h*d[i];
  28.         for(i=0; i<n; i++)
  29. { 
  30. y[i]=(y[i]+d[i])/2.0;
  31.             z(i,j)=y[i];
  32. }
  33. }
  34. }
  36. //变步长积分欧拉法
  37. template <class _Ty>
  38. void ODE_EulerVariationalStep(_Ty t, _Ty h, valarray<_Ty>& y, _Ty eps)
  39. {
  40. int j, m(1);
  41.     _Ty p(1.0+eps), x, q, hh=h;
  43. int n = y.size(); //微分方程组中方程的个数，也是未知函数的个数
  44.     valarray<_Ty> a(n), b(n), c(n), d(n);
  46.     for(int i=0; i<n; i++) a[i] = y[i];
  48.     while(p>eps || FloatEqual(p,eps))
  49. { 
  50. for(i=0; i<n; i++)
  51. {
  52. b[i]=y[i];
  53. y[i]=a[i];
  54. }
  55.         for(j=0; j<m; j++)
  56. {
  57. for(i=0; i<n; i++) c[i]=y[i];
  58.             x=t+j*hh;
  59.             FunctionValueEVS(x,y,d);
  60.             for(i=0; i<n; i++)
  61. y[i]=c[i]+hh*d[i];
  62.             x=t+(j+1)*hh;
  63.             FunctionValueEVS(x,y,d);
  64.             for(i=0; i<n; i++)
  65. d[i]=c[i]+hh*d[i];
  66.             for(i=0; i<n; i++)
  67. y[i]=(y[i]+d[i])/2.0;
  68. }
  69.         p = 0.0;
  70.         for(i=0; i<n; i++)
  71. { 
  72. q=Abs(y[i]-b[i]);
  73.             if(q>p) p=q;
  74.         }
  75.         hh /= 2.0;
  76. m = m + m;
  77.     }
  78. }
  80. //全区间定步长积分维梯法
  81. template <class _Ty>
  82. void ODE_WittyContentStep(_Ty t, valarray<_Ty>& y,  
  83. _Ty h, int k, matrix<_Ty>& z)
  84. {
  85.     int n = y.size(); //微分方程组中方程的个数，也是未知函数的个数
  86. valarray<_Ty> a(n), d(n);
  88.     for(int i=0; i<n; i++) z(i,0)=y[i];
  89.     FunctionValueWCS(t,y,d);
  90.     for(int j=1; j<k; j++)
  91.     { 
  92. for(i=0; i<n; i++)
  93. a[i]=z(i,j-1)+h*d[i]/2.0;
  94.         _Ty x=t+(j-0.5)*h;
  95.         FunctionValueWCS(x,a,y);
  96.         for(i=0; i<n; i++)
  97.         { 
  98. d[i]=2.0*y[i]-d[i];
  99.             z(i,j)=z(i,j-1)+h*y[i];
  100.         }
  101.     }
  102. }
  104. //全区间定步长积分龙格-库塔法
  105. template <class _Ty>
  106. void ODE_RungeKuttaContentStep(_Ty t, valarray<_Ty>& y,  
  107. _Ty h, int k, matrix<_Ty>& z)
  108. {
  109.     _Ty a[4], tt;
  111. int n = y.size(); //微分方程组中方程的个数，也是未知函数的个数
  112. valarray<_Ty> b(n), d(n);
  114.     a[0] = h / 2.0;
  115. a[1] = a[0];
  116.     a[2] = a[3] = h;
  118.     for(int i=0; i<n; i++) z(i,0)=y[i];
  120.     for(int l=1; l<k; l++)
  121.     { 
  122. FunctionValueRKCS(t, y, d);
  123.         for(i=0; i<n; i++) b[i]=y[i];
  124.         for(int j=0; j<=2; j++)
  125.         {
  126. for(i=0; i<n; i++)
  127.             {
  128. y[i]=z(i,l-1)+a[j]*d[i];
  129.                 b[i]=b[i]+a[j+1]*d[i]/3.0;
  130.             }
  131.             tt=t+a[j];
  132.             FunctionValueRKCS(tt, y, d);
  133.         }
  134.         for(i=0; i<n; i++)
  135.           y[i]=b[i]+h*d[i]/6.0;
  136.         for(i=0; i<n; i++)
  137.           z(i,l)=y[i];
  138.         t=t+h;
  139.     }
  140. }
  142. //变步长积分龙格-库塔法
  143. template <class _Ty>
  144. void ODE_RungeKuttaVariationalStep(_Ty t, _Ty h, 
  145. valarray<_Ty>& y, _Ty eps)
  146. {
  147. int m(1), j, k;
  148.     _Ty hh(h), dt,x(t),tt,q,a[4];
  150. int n = y.size(); //微分方程组中方程的个数，也是未知函数的个数
  151. valarray<_Ty> g(n), b(n), c(n), d(n), e(n);
  152.     _Ty p=1.0+eps;
  154.     for(int i=0; i<n; i++) c[i]=y[i];
  155.     while(p>eps || FloatEqual(p,eps))
  156.     {
  157. a[0]=hh/2.0; 
  158. a[1]=a[0];
  159. a[2] = a[3] = hh;
  160.         for(i=0; i<n; i++)
  161.         {
  162. g[i]=y[i];
  163. y[i]=c[i];
  164. }
  165.         dt=h/m; 
  166. t=x;
  167.         for(j=0; j<m; j++)
  168.         {
  169. FunctionValueRKVS(t, y, d);
  170.             for(i=0; i<n; i++) 
  171.             {
  172. b[i]=y[i]; 
  173. e[i]=y[i];
  174. }
  175.             for(k=0; k<=2; k++)
  176.             {
  177. for(i=0; i<n; i++)
  178. {
  179. y[i]=e[i]+a[k]*d[i];
  180.                     b[i]=b[i]+a[k+1]*d[i]/3.0;
  181.                 }
  182.                 tt=t+a[k];
  183.                 FunctionValueRKVS(tt, y, d);
  184.             }
  185.             for(i=0; i<n; i++)
  186.               y[i]=b[i]+hh*d[i]/6.0;
  187.             t += dt;
  188.         }
  189.         p = 0.0;
  190.         for(i=0; i<n; i++)
  191.         {
  192. q = Abs(y[i]-g[i]);
  193.             if(q>p) p = q;
  194.         }
  195.         hh /= 2.0;
  196. m = m + m;
  197.     }
  198. }
  200. //变步长积分基尔法
  201. template <class _Ty>
  202. void ODE_GillVariationalStep(_Ty t, _Ty h, valarray<_Ty>& y,
  203. _Ty eps, valarray<_Ty>& q)
  204. {
  205.     int i,k,m(1),ii;
  206.     _Ty x(t),hh(h),r,s,t0,dt,qq;
  207.     _Ty a[4]={0.5, 0.29289321881, 1.7071067812, 0.166666667};
  208.     _Ty b[4]={2.0,1.0,1.0,2.0};
  209.     _Ty c[4],e[4]={0.5,0.5,1.0,1.0};
  211. int n = y.size(); //微分方程组中方程的个数，也是未知函数的个数
  212. valarray<_Ty> g(n), u(n), v(n), d(n);
  214.     for(i=0; i<3; i++) c[i]=a[i];
  215.     c[3]=0.5;
  216.     _Ty p=1.0+eps;
  218.     for(int j=0; j<n; j++) u[j]=y[j];
  219.     while(p>eps || FloatEqual(p,eps))
  220.     {
  221. for(j=0; j<n; j++)
  222. {
  223. v[j]=y[j];
  224. y[j]=u[j];
  225. g[j]=q[j];
  226. }
  227.         dt=h/m; 
  228. t=x;
  229.         for(k=0; k<m; k++)
  230.         { 
  231. FunctionValueGVS(t,y,d);
  232.             for(ii=0; ii<=3; ii++)
  233.             {
  234. for(j=0; j<n; j++)
  235.                   d[j]=d[j]*hh;
  236.                 for(j=0; j<n; j++)
  237.                 {
  238. r=(a[ii]*(d[j]-b[ii]*g[j])+y[j])-y[j];
  239.                     y[j]=y[j]+r;
  240.                     s=g[j]+3.0*r;
  241.                     g[j]=s-c[ii]*d[j];
  242.                 }
  243.                 t0=t+e[ii]*hh;
  244.                 FunctionValueGVS(t0,y,d);
  245.             }
  246.             t=t+dt;
  247.         }
  248.         p=0.0;
  249.         for(j=0; j<n; j++)
  250.         {
  251. qq=Abs(y[j]-v[j]);
  252.             if(qq>p) p=qq;
  253.         }
  254.         hh=hh/2.0; 
  255. m=m+m;
  256.     }
  257.     for(j=0; j<n; j++) q[j]=g[j];
  258. }
  260. //全区间变步长积分基尔法
  261. template <class _Ty>
  262. void ODE_GillVariationalStepInterval(_Ty t, _Ty h, valarray<_Ty>& y,
  263. _Ty eps, int k, matrix<_Ty>& z)
  264. {
  265.     int j,m,kk,ii;
  266.     _Ty aa(t),hh,x,p,dt,r,s,t0,qq;
  267.     _Ty a[4]={0.5, 0.29289321881, 1.7071067812, 0.166666667};
  268.     _Ty b[4]={2.0,1.0,1.0,2.0};
  269.     _Ty c[4],e[4]={0.5,0.5,1.0,1.0};
  271. int n = y.size(); //微分方程组中方程的个数，也是未知函数的个数
  272. valarray<_Ty> g(n), q(n), u(n), v(n), d(n);
  273.     
  274.     for(int i=0; i<3; i++) c[i]=a[i];
  275.     c[3]=0.5;
  276.     
  277.     for(i=0; i<n; i++) 
  278.     {
  279. z(i,0) = y[i];
  280. q[i] = 0.0;
  281. }
  282.     for(i=2; i<=k; i++)
  283.     {
  284. x=aa+(i-2)*h;
  285. m=1;
  286. hh=h;
  287.         p=1.0+eps;
  288.         for(j=0; j<n; j++) u[j]=y[j];
  289.         while(p>eps || FloatEqual(p,eps))
  290.         {
  291. for(j=0; j<n; j++)
  292. { 
  293. v[j]=y[j];
  294. y[j]=u[j];
  295. g[j]=q[j];
  296. }
  297.             dt=h/m; 
  298. t=x;
  299.             for(kk=0; kk<m; kk++)
  300.             {
  301. FunctionValueGVSI(t,y,d);
  302.                 for(ii=0; ii<4; ii++)
  303.                 {
  304. for(j=0; j<n; j++)
  305.                       d[j]=d[j]*hh;
  306.                     for(j=0; j<n; j++)
  307.                     {
  308. r=(a[ii]*(d[j]-b[ii]*g[j])+y[j])-y[j];
  309.                         y[j]=y[j]+r;
  310.                         s=g[j]+3.0*r;
  311.                         g[j]=s-c[ii]*d[j];
  312.                     }
  313.                     t0=t+e[ii]*hh;
  314.                     FunctionValueGVSI(t0,y,d);
  315.                 }
  316.                 t += dt;
  317.             }
  318.             p=0.0;
  319.             for(j=0; j<n; j++)
  320.             { 
  321. qq=Abs(y[j]-v[j]);
  322.                 if(qq>p) p=qq;
  323.             }
  324.             hh /= 2.0; 
  325. m=m+m;
  326.         }
  327.         for(j=0; j<n; j++)
  328.         { 
  329. q[j]=g[j];
  330. z(j,i-1)=y[j];
  331. }
  332.     }
  333. }
  335. //变步长积分默森法
  336. template <class _Ty>
  337. void ODE_MersonVariationalStep(_Ty t, _Ty h, int n, valarray<_Ty>& y,
  338. _Ty eps, int k, matrix<_Ty>& z)
  339. {
  340.     int j,m,nn;
  341.     _Ty aa(t),bb,x,hh,p,dt,t0,qq;
  342.     valarray<_Ty> a(n), b(n), c(n), u(n), v(n), d(n);
  343.         
  344.     for(int i=0; i<n; i++) z(i,0)=y[i];
  345.     for(i=1; i<k; i++)
  346.     { 
  347. x=aa+(i-1)*h; 
  348. nn=1; 
  349. hh=h;
  350.         for(j=0; j<n; j++) u[j]=y[j];
  351.         p=1.0+eps;
  352.         while(p>eps || FloatEqual(p,eps))
  353.         {
  354. for(j=0; j<n; j++)
  355. {
  356. v[j]=y[j]; 
  357. y[j]=u[j];
  358. }
  359.             dt=h/nn; t=x;
  360.             for(m=0; m<nn; m++)
  361. { 
  362. FunctionValueMVS(t,y,d);
  363.                 for(j=0; j<n; j++)
  364.                 {
  365. a[j]=d[j];
  366. y[j]=y[j]+hh*d[j]/3.0;
  367. }
  368.                 t0=t+hh/3.0;
  369.                 FunctionValueMVS(t0,y,d);
  370.                 for(j=0; j<n; j++)
  371.                 {
  372. b[j]=d[j];
  373. y[j]=y[j]+hh*(d[j]-a[j])/6.0;
  374. }
  375.                 FunctionValueMVS(t0,y,d);
  376.                 for(j=0; j<n; j++)
  377.                 {
  378. b[j]=d[j];
  379.                     bb=(d[j]-4.0*(b[j]+a[j]/4.0)/9.0)/8.0;
  380.                     y[j]=y[j]+3.0*hh*bb;
  381.                 }
  382.                 t0=t+hh/2.0;
  383.                 FunctionValueMVS(t0,y,d);
  384.                 for(j=0; j<n; j++)
  385.                 { c[j]=d[j];
  386.                     qq=d[j]-15.0*(b[j]-a[j]/5.0)/16.0;
  387.                     y[j]=y[j]+2.0*hh*qq;
  388.                 }
  389.                 t0=t+hh;
  390.                 FunctionValueMVS(t0,y,d);
  391.                 for(j=0; j<n; j++)
  392.                 { 
  393. qq=c[j]-9.0*(b[j]-2.0*a[j]/9.0)/8.0;
  394.                     qq=d[j]-8.0*qq;
  395.                     y[j]=y[j]+hh*qq/6.0;
  396.                 }
  397.                 t += dt;
  398.             }
  399.             p=0.0;
  400.             for(j=0; j<n; j++)
  401.             { qq=Abs(y[j]-v[j]);
  402.                 if(qq>p) p=qq;
  403.             }
  404.             hh /= 2.0;
  405. nn=nn+nn;
  406.         }
  407.         for(j=0; j<n; j++) z(j,i)=y[j];
  408.     }
  409. }
  411. //积分一步连分式法
  412. template <class _Ty>
  413. void ODE_Fraction(_Ty t, _Ty h, valarray<_Ty>& y, _Ty eps)
  414. {
  415. int i,j,k(1),m,nn(1),it(1);
  417. int n = y.size(); //微分方程组中方程的个数，也是未知函数的个数
  418.     _Ty x(t),hh(h),dd,q,p,g[10];
  419. valarray<_Ty> e(n), b(10*n), w(n), u(n), v(n), d(n);
  421.     for(j=0; j<n; j++) v[j]=y[j];
  422. g[0]=hh;
  423.     Fraction1(x,hh,y,w,d,e);
  424.     for(j=0; j<n; j++)
  425.     { 
  426. b[j]=y[j];
  427. u[j]=y[j];
  428. }
  430.     while(it==1)
  431.     { 
  432. nn=nn+nn; 
  433. hh=hh/2.0;
  434. it=0;
  435.         g[k]=hh;
  436.         for(j=0; j<n; j++) y[j]=v[j];
  437.         t=x;
  438.         for(j=0; j<nn; j++)
  439.         { 
  440. Fraction1(t,hh,y,w,d,e);
  441.             t=t+hh;
  442.         }
  443.         for(j=0; j<n; j++)
  444.         { 
  445. dd=y[j]; 
  446. m=0;
  447.             for(i=0; i<k; i++)
  448.               if(m==0)
  449.               {
  450.   q=dd-b[i*n+j];
  451.                   if(FloatEqual(q,0.0)) m=1;
  452.                   else dd=(g[k]-g[i])/q;
  453.               }
  454.             b[k*n+j]=dd;
  455.             if(m!=0) b[k*n+j]=1.0e+35;
  456.         }
  457.         for(j=0; j<n; j++)
  458.         {
  459. dd=0.0;
  460.             for(i=k-1; i>=0; i--)
  461. dd=-g[i]/(b[(i+1)*n+j]+dd);
  462.             y[j]=dd+b[j];
  463.         }
  464.         p=0.0;
  465.         for(j=0; j<n; j++)
  466.         { 
  467. q=Abs(y[j]-u[j]);
  468.             if(q>p) p=q;
  469.         }
  470.         if((p>eps || FloatEqual(p,eps))&&(k<7))
  471.         { 
  472. for(j=0; j<n; j++) u[j]=y[j];
  473.             k=k+1; 
  474. it=1;
  475.         }
  476.     }
  477. }
  479. //积分一步连分式法ODE_Fraction()辅助函数
  480. template <class _Ty>
  481. void Fraction1(_Ty t, _Ty h, valarray<_Ty>& y, valarray<_Ty>& b, 
  482. valarray<_Ty>& d, valarray<_Ty>& e)
  483. {
  484.     _Ty a[4],tt;
  485. int n = y.size(); //微分方程组中方程的个数，也是未知函数的个数
  486.     
  487. a[0]=h/2.0;
  488. a[1]=a[0];
  489. a[2]=h; 
  490. a[3]=h;
  491.     FunctionValueF(t,y,d);
  492.     for(int i=0; i<n; i++)
  493. {
  494. b[i]=y[i];
  495. e[i]=y[i];
  496. }
  497.     for(int k=0; k<=2; k++)
  498.     {
  499. for(i=0; i<n; i++)
  500.         {
  501. y[i]=e[i]+a[k]*d[i];
  502.             b[i]=b[i]+a[k+1]*d[i]/3.0;
  503.         }
  504.         tt=t+a[k];
  505.         FunctionValueF(tt,y,d);
  506.     }
  507.     for(i=0; i<n; i++)
  508.       y[i]=b[i]+h*d[i]/6.0;
  509. }
  512. //全区间积分连分式法
  513. template <class _Ty>
  514. void ODE_FractionInterval(_Ty t, _Ty h, valarray<_Ty>& y,
  515. _Ty eps, int k, matrix<_Ty>& z)
  516. {
  517. int j,kk,l,m,nn,it;
  518.     _Ty x,hh,dd,tt,q,p,g[10];
  520. int n = y.size(); //微分方程组中方程的个数，也是未知函数的个数
  521. valarray<_Ty> e(n), b(10*n), w(n), u(n), v(n), d(n);
  523.     for(int i=0; i<n; i++) z(i,0)=y[i];
  524.     for(i=1; i<k; i++)
  525.     { 
  526. for(j=0; j<n; j++) v[j]=y[j];
  527.         x=t+(i-1)*h; 
  528. nn=1; 
  529. hh=h; 
  530. g[0]=hh;
  531.         Fraction2(x,hh,y,w,d,e);
  532.         for(j=0; j<n; j++)
  533.         {
  534. b[j]=y[j];
  535. u[j]=y[j];
  536. }
  537.         kk=1;
  538. it=1;
  539.         while(it==1)
  540.         {
  541. nn=nn+nn; 
  542. hh=hh/2.0;
  543. it=0;
  544.             g[kk]=hh;
  545.             for(j=0; j<n; j++) y[j]=v[j];
  546.             tt=x;
  547.             for(j=0; j<nn; j++)
  548.             {
  549. Fraction2(t,hh,y,w,d,e);
  550.                 tt=tt+hh;
  551.             }
  552.             for(j=0; j<n; j++)
  553.             { 
  554. dd=y[j]; 
  555. l=0;
  556.                 for(m=0; m<kk; m++)
  557.                   if(l==0)
  558.                   {
  559.   q=dd-b[m*n+j];
  560.                       if(Abs(q)+1.0==1.0) l=1;
  561.                       else dd=(g[kk]-g[m])/q;
  562.                   }
  563.                 b[kk*n+j]=dd;
  564.                 if(l!=0) b[kk*n+j]=1.0e+35;
  565.             }
  566.             for(j=0; j<n; j++)
  567.             { 
  568. dd=0.0;
  569.                 for(m=kk-1; m>=0; m--)
  570. dd=-g[m]/(b[(m+1)*n+j]+dd);
  571.                 y[j]=dd+b[j];
  572.             }
  573.             p=0.0;
  574.             for(j=0; j<n; j++)
  575.             { 
  576. q=Abs(y[j]-u[j]);
  577.                 if(q>p) p=q;
  578.             }
  579.             if((p>eps || FloatEqual(p,eps))&&(kk<7))
  580.             { for(j=0; j<n; j++) u[j]=y[j];
  581.                 kk=kk+1; it=1;
  582.             }
  583.         }
  584.         for(j=0; j<n; j++)
  585.           z(j,i)=y[j];
  586.     }
  587. }
  589. //全区间积分连分式法辅助函数
  590. template <class _Ty>
  591. void Fraction2(_Ty t, _Ty h, valarray<_Ty>& y, valarray<_Ty>& b, 
  592. valarray<_Ty>& d, valarray<_Ty>& e)
  593. {
  594.     int i,k;
  595.     _Ty a[4],tt;
  597. int n = y.size(); //微分方程组中方程的个数，也是未知函数的个数
  598.     a[0]=h/2.0; 
  599. a[1]=a[0];
  600. a[2]=h;
  601. a[3]=h;
  602.     FunctionValueFI(t,y,d);
  603.     for(i=0; i<n; i++) 
  604. {
  605. b[i]=y[i]; 
  606. e[i]=y[i];
  607. }
  608.     for(k=0; k<3; k++)
  609.     { 
  610. for(i=0; i<n; i++)
  611.         {
  612. y[i]=e[i]+a[k]*d[i];
  613.             b[i]=b[i]+a[k+1]*d[i]/3.0;
  614.         }
  615.         tt=t+a[k];
  616.         FunctionValueFI(tt,y,d);
  617.     }
  618.     for(i=0; i<n; i++)
  619.       y[i]=b[i]+h*d[i]/6.0;
  620. }
  622. //全区间积分双边法
  623. template <class _Ty>
  624. void ODE_BothSidesInterval(_Ty t, _Ty h, valarray<_Ty>& y,
  625.     _Ty eps, int k, matrix<_Ty>& z)
  626. {
  627.     _Ty a(t), qq;
  629. int n = y.size(); //微分方程组中方程的个数，也是未知函数的个数
  630. valarray<_Ty> d(n), w(n), u(n), v(n), p(n);
  632.     for(int i=0; i<n; i++)
  633. {
  634. p[i]=0.0;
  635. z(i,0)=y[i];
  636. }
  638.     FunctionValueBSI(t,y,d);
  639.     for(int j=0; j<n; j++) u[j]=d[j];
  640.     BothSidesInterval(t,h,y,eps); //辅助函数
  641.     t=a+h;
  642.     FunctionValueBSI(t,y,d);
  643.     for(j=0; j<n; j++)
  644.     { 
  645. z(j,1)=y[j]; v[j]=d[j];}
  646.     for(j=0; j<n; j++)
  647.     { 
  648. p[j]=-4.0*z(j,1)+5.0*z(j,0)+2.0*h*(2.0*v[j]+u[j]);
  649.         y[j]=p[j];
  650.     }
  651.     t=a+2.0*h;
  652.     FunctionValueBSI(t,y,d);
  653.     for(j=0; j<n; j++)
  654.     { 
  655. qq=2.0*h*(d[j]-2.0*v[j]-2.0*u[j])/3.0;
  656.         qq=qq+4.0*z(j,1)-3.0*z(j,0);
  657.         z(j,2)=(p[j]+qq)/2.0;
  658.         y[j]=z(j,2);
  659.     }
  660.     for(i=3; i<k; i++)
  661.     { 
  662. t=a+(i-1)*h;
  663.         FunctionValueBSI(t,y,d);
  664.         for(j=0; j<n; j++)
  665.         {
  666. u[j]=v[j];
  667. v[j]=d[j];
  668. }
  669.         for(j=0; j<n; j++)
  670.         { 
  671. qq=-4.0*z(j,i-1)+5.0*z(j,i-2);
  672.             p[j]=qq+2.0*h*(2.0*v[j]+u[j]);
  673.             y[j]=p[j];
  674.         }
  675.         t=t+h;
  676.         FunctionValueBSI(t,y,d);
  677.         for(j=0; j<n; j++)
  678.         { 
  679. qq=2.0*h*(d[j]-2.0*v[j]-2.0*u[j])/3.0;
  680.             qq=qq+4.0*z(j,i-1)-3.0*z(j,i-2);
  681.             y[j]=(p[j]+qq)/2.0;
  682.             z(j,i)=y[j];
  683.         }
  684.     }
  685. }
  687. //全区间积分双边法辅助函数
  688. template <class _Ty>
  689. void BothSidesInterval(_Ty t, _Ty h, valarray<_Ty>& y, _Ty eps)
  690. { 
  691.     int m(1), j, k;
  692.     _Ty hh(h), dt,x(t), tt, q, a[4];
  694. int n = y.size(); //微分方程组中方程的个数，也是未知函数的个数
  695. valarray<_Ty> d(n), b(n), c(n), e(n), g(n);
  696.     
  697.     _Ty p=1.0+eps;
  698.     for(int i=0; i<n; i++) c[i]=y[i];
  699.     while(p>eps || FloatEqual(p,eps))
  700.     { 
  701. a[0]=hh/2.0;
  702. a[1]=a[0]; 
  703. a[2]=hh;
  704. a[3]=hh;
  705.         for(i=0; i<n; i++)
  706.         {
  707. g[i]=y[i]; 
  708. y[i]=c[i];
  709. }
  710.         dt=h/m; 
  711. t=x;
  712.         for(j=0; j<m; j++)
  713.         {
  714. FunctionValueBSI(t,y,d);
  715.             for(i=0; i<n; i++) 
  716.             {
  717. b[i]=y[i];
  718. e[i]=y[i];
  719. }
  720.             for(k=0; k<3; k++)
  721.             { 
  722. for(i=0; i<n; i++)
  723.                 {
  724. y[i]=e[i]+a[k]*d[i];
  725.                     b[i]=b[i]+a[k+1]*d[i]/3.0;
  726.                 }
  727.                 tt=t+a[k];
  728.                 FunctionValueBSI(tt,y,d);
  729.             }
  730.             for(i=0; i<n; i++)
  731.               y[i]=b[i]+hh*d[i]/6.0;
  732.             t=t+dt;
  733.         }
  734.         p=0.0;
  735.         for(i=0; i<n; i++)
  736.         { 
  737. q=Abs(y[i]-g[i]);
  738.             if(q>p) p=q;
  739.         }
  740.         hh=hh/2.0;
  741. m=m+m;
  742.     }
  743. }
  745. //全区间积分阿当姆斯预报-校正法
  746. template <class _Ty>
  747. void ODE_AdamsInterval(_Ty t, _Ty h, valarray<_Ty>& y,
  748. _Ty eps, int k, matrix<_Ty>& z)
  749. {
  750.     int j,m;
  751.     _Ty a(t), q;
  753. int n = y.size(); //微分方程组中方程的个数，也是未知函数的个数
  754. valarray<_Ty> b(4*n), e(n), s(n), g(n), d(n);
  755.     
  756.     for(int i=0; i<n; i++) z(i,0)=y[i];
  757.     FunctionValueAI(t,y,d);
  758.     for(i=0; i<n; i++) b[i]=d[i];
  759.     for(i=1; i<4; i++)
  760.       if(i<k)
  761.       { 
  762.   t=a+i*h;
  763.           AdamsInterval(t,h,y,eps); //辅助函数
  764.           for(j=0; j<n; j++) z(j,i)=y[j];
  765.           FunctionValueAI(t,y,d);
  766.           for(j=0; j<n; j++) b[i*n+j]=d[j];
  767.       }
  768.     for(i=4; i<k; i++)
  769.     { 
  770. for(j=0; j<n; j++)
  771.         { 
  772. q=55.0*b[3*n+j]-59.0*b[2*n+j];
  773.             q=q+37.0*b[n+j]-9.0*b[j];
  774.             y[j]=z(j,i-1)+h*q/24.0;
  775.             b[j]=b[n+j];
  776.             b[n+j]=b[n+n+j];
  777.             b[n+n+j]=b[n+n+n+j];
  778.         }
  779.         t=a+i*h;
  780.         FunctionValueAI(t,y,d);
  781.         for(m=0; m<n; m++) b[n+n+n+m]=d[m];
  782.         for(j=0; j<n; j++)
  783.         { 
  784. q=9.0*b[3*n+j]+19.0*b[n+n+j]-5.0*b[n+j]+b[j];
  785.             y[j]=z(j,i-1)+h*q/24.0;
  786.             z(j,i)=y[j];
  787.         }
  788.         FunctionValueAI(t,y,d);
  789.         for(m=0; m<n; m++) b[3*n+m]=d[m];
  790.     }
  791. }
  793. //全区间积分阿当姆斯预报-校正法辅助函数
  794. template <class _Ty>
  795. void AdamsInterval(_Ty t, _Ty h, valarray<_Ty>& y, _Ty eps)
  796. {
  797.     int m(1), j, k;
  799. int n = y.size(); //微分方程组中方程的个数，也是未知函数的个数
  800.     _Ty hh(h), dt, x(t), tt, q, a[4];
  801. valarray<_Ty> g(n), b(n), c(n), d(n), e(n);
  802.     
  803.     _Ty p=1.0+eps;
  804.     for(int i=0; i<n; i++) c[i]=y[i];
  805.     while(p>eps || FloatEqual(p,eps))
  806.     { 
  807. a[0]=hh/2.0;
  808. a[1]=a[0]; 
  809. a[2]=hh;
  810. a[3]=hh;
  811.         for(i=0; i<n; i++)
  812.         {
  813. g[i]=y[i]; 
  814. y[i]=c[i];
  815. }
  816.         dt=h/m; t=x;
  817.         for(j=0; j<m; j++)
  818.         { 
  819. FunctionValueAI(t,y,d);
  820.             for(i=0; i<n; i++) 
  821.             {
  822. b[i]=y[i];
  823. e[i]=y[i];
  824. }
  825.             for(k=0; k<3; k++)
  826.             { 
  827. for(i=0; i<n; i++)
  828.                 {
  829. y[i]=e[i]+a[k]*d[i];
  830.                     b[i]=b[i]+a[k+1]*d[i]/3.0;
  831.                 }
  832.                 tt=t+a[k];
  833.                 FunctionValueAI(tt,y,d);
  834.             }
  835.             for(i=0; i<n; i++)
  836.               y[i]=b[i]+hh*d[i]/6.0;
  837.             t=t+dt;
  838.         }
  839.         p=0.0;
  840.         for(i=0; i<n; i++)
  841.         {
  842. q=Abs(y[i]-g[i]);
  843.             if(q>p) p=q;
  844.         }
  845.         hh=hh/2.0;
  846. m=m+m;
  847.     }
  848. }
  850. //全区间积分哈明法
  851. template <class _Ty>
  852. void ODE_HammingInterval(_Ty t, _Ty h, valarray<_Ty>& y,
  853. _Ty eps, int k, matrix<_Ty>& z)
  854. {
  855.     int j,m;
  856.     _Ty a(t), q;
  858. int n = y.size(); //微分方程组中方程的个数，也是未知函数的个数
  859. valarray<_Ty> d(n), w(n), u(n), v(n), b(4*n), g(n);
  860.     
  861.     for(int i=0; i<n; i++) z(i,0)=y[i];
  862.     FunctionValueHI(t,y,d);
  863.     for(i=0; i<n; i++) b[i]=d[i];
  864.     for(i=1; i<=3; i++)
  865.       if(i<k)
  866.       { 
  867.   t=a+i*h;
  868.           HammingInterval(t,h,y,eps);
  869.           for(m=0; m<n; m++) z(m,i)=y[m];
  870.           FunctionValueHI(t,y,d);
  871.           for(m=0; m<n; m++) b[i*n+m]=d[m];
  872.       }
  873.     for(i=0; i<n; i++) u[i]=0.0;
  874.     for(i=4; i<k; i++)
  875.     { 
  876. for(j=0; j<n; j++)
  877.         {
  878. q=2.0*b[3*n+j]-b[n+n+j]+2.0*b[n+j];
  879.             y[j]=z(j,i-4)+4.0*h*q/3.0;
  880.         }
  881.         for(j=0; j<n; j++)
  882. y[j]=y[j]+112.0*u[j]/121.0;
  883.         t=a+i*h;
  884.         FunctionValueHI(t,y,d);
  885.         for(j=0; j<n; j++)
  886.         { 
  887. q=9.0*z(j,i-1)-z(j,i-3);
  888.             q=(q+3.0*h*(d[j]+2.0*b[3*n+j]-b[n+n+j]))/8.0;
  889.             u[j]=q-y[j];
  890.             z(j,i)=q-9.0*u[j]/121.0;
  891.             y[j]=z(j,i);
  892.             b[n+j]=b[n+n+j];
  893.             b[n+n+j]=b[n+n+n+j];
  894.         }
  895.         FunctionValueHI(t,y,d);
  896.         for(m=0; m<n; m++) b[3*n+m]=d[m];
  897.     }
  898. }
  900. //全区间积分哈明法辅助函数
  901. template <class _Ty>
  902. void HammingInterval(_Ty t, _Ty h, valarray<_Ty>& y, _Ty eps)
  903. { 
  904.     int m(1), j, k;
  905.     _Ty hh(h), dt, x(t), tt, q, a[4];
  907. int n = y.size(); //微分方程组中方程的个数，也是未知函数的个数
  908. valarray<_Ty> g(n), b(n), c(n), d(n), e(n);
  909.     
  910. _Ty p=1.0+eps;
  911.     for(int i=0; i<n; i++) c[i]=y[i];
  912.     while(p>eps || FloatEqual(p,eps))
  913.     { 
  914. a[0]=hh/2.0; 
  915. a[1]=a[0]; 
  916. a[2]=hh;
  917. a[3]=hh;
  918.         for(i=0; i<n; i++)
  919.         {
  920. g[i]=y[i]; 
  921. y[i]=c[i];
  922. }
  923.         dt=h/m;
  924. t=x;
  925.         for(j=0; j<m; j++)
  926.         {
  927. FunctionValueHI(t,y,d);
  928.             for(i=0; i<n; i++) 
  929.             {
  930. b[i]=y[i];
  931. e[i]=y[i];
  932. }
  933.             for(k=0; k<3; k++)
  934.             { 
  935. for(i=0; i<n; i++)
  936.                 {
  937. y[i]=e[i]+a[k]*d[i];
  938.                     b[i]=b[i]+a[k+1]*d[i]/3.0;
  939.                 }
  940.                 tt=t+a[k];
  941.                 FunctionValueHI(t,y,d);
  942.             }
  943.             for(i=0; i<n; i++)
  944. y[i]=b[i]+hh*d[i]/6.0;
  945.             t=t+dt;
  946.         }
  947.         p=0.0;
  948.         for(i=0; i<n; i++)
  949.         { 
  950. q=Abs(y[i]-g[i]);
  951.             if(q>p) p=q;
  952.         }
  953.         hh=hh/2.0; 
  954. m=m+m;
  955.     }
  956. }
  958. //积分一步特雷纳法
  959. template <class _Ty>
  960. void ODE_Treanor(_Ty t, _Ty h, valarray<_Ty>& y)
  961. {
  962. //extern void gtnr1f();
  963.     _Ty s,aa,bb,dd,g,dy,dy1;
  965. int n = y.size(); //微分方程组中方程的个数，也是未知函数的个数
  966. valarray<_Ty> w(4*n), q(4*n), r(4*n), d(n), p(n);
  967.     
  968.     for(int j=0; j<n; j++) w[j]=y[j];
  969.     FunctionValueT(t,y,d);
  970.     for(j=0; j<n; j++)
  971.     {
  972. q[j]=d[j]; 
  973. y[j]=w[j]+h*d[j]/2.0;
  974.         w[n+j]=y[j];
  975.     }
  976.     s=t+h/2.0;
  977.     FunctionValueT(t,y,d);
  978.     for(j=0; j<n; j++)
  979.     { 
  980. q[n+j]=d[j];
  981.         y[j]=w[j]+h*d[j]/2.0;
  982.         w[n+n+j]=y[j];
  983.     }
  984.     FunctionValueT(t,y,d);
  985.     for(j=0; j<n; j++) q[n+n+j]=d[j];
  987.     for(j=0; j<n; j++)
  988.     { 
  989. aa=q[n+n+j]-q[n+j];
  990.         bb=w[n+n+j]-w[n+j];
  991.         if(-aa*bb*h>0.0)
  992.         {
  993. p[j]=-aa/bb; dd=-p[j]*h;
  994.             r[j]=exp(dd);
  995.             r[n+j]=(r[j]-1.0)/dd;
  996.             r[n+n+j]=(r[n+j]-1.0)/dd;
  997.             r[3*n+j]=(r[n+n+j]-1.0)/dd;
  998.         }
  999.         else p[j]=0.0;
  1000.         if(p[j]<0.0||FloatEqual(p[j],0.0)) g=q[n+n+j];
  1001.         else
  1002.         { 
  1003. g=2.0*(q[n+n+j]-q[j])*r[n+n+j];
  1004.             g=g+(q[j]-q[n+j])*r[n+j]+q[n+j];
  1005.         }
  1006.         w[3*n+j]=w[j]+g*h;
  1007.         y[j]=w[3*n+j];
  1008.     }
  1009.     s=t+h;
  1010.     FunctionValueT(t,y,d);
  1011.     for(j=0; j<n; j++) q[3*n+j]=d[j];
  1013.     for(j=0; j<n; j++)
  1014.     { 
  1015. if(p[j]<0.0||FloatEqual(p[j],0.0))
  1016.         { 
  1017. dy=q[j]+2.0*(q[n+j]+q[n+n+j]);
  1018.             dy=(dy+q[n+n+n+j])*h/6.0;
  1019.         }
  1020.         else
  1021.         {
  1022. dy=-3.0*(q[j]+p[j]*w[j])+2.0*(q[n+j]+p[j]*w[n+j]);
  1023.             dy=dy+2.0*(q[n+n+j]+p[j]*w[n+n+j]);
  1024.             dy=dy-(q[n+n+n+j]+p[j]*w[n+n+n+j]);
  1025.             dy=dy*r[n+n+j]+q[j]*r[n+j];
  1026.             dy1=q[j]-q[n+j]-q[n+n+j]+q[n+n+n+j];
  1027.             dy1=dy1+(w[j]-w[n+j]-w[n+n+j]+w[n+n+n+j])*p[j];
  1028.             dy=(dy+4.0*dy1*r[n+n+n+j])*h;
  1029.         }
  1030.         y[j]=w[j]+dy;
  1031.     }
  1032. }
  1034. //全区间积分特雷纳法
  1035. template <class _Ty>
  1036. void ODE_TreanorInterval(_Ty t, _Ty h, valarray<_Ty>& y,
  1037. int k, matrix<_Ty>& z)
  1038. {
  1039.     _Ty a(t),s,aa,bb,dd,g,dy,dy1;
  1041. int n = y.size(); //微分方程组中方程的个数，也是未知函数的个数
  1042. valarray<_Ty> w(4*n), q(4*n), r(4*n), d(n), p(n);
  1043.     
  1044.     for(int j=0; j<n; j++) z(j,0)=y[j];
  1046.     for(int i=1; i<k; i++)
  1047.     {
  1048. t=a+(i-1)*h;
  1049.         for(j=0; j<n; j++) w[j]=y[j];
  1050.         FunctionValueTI(s,y,d);
  1051.         for(j=0; j<n; j++)
  1052.         { 
  1053. q[j]=d[j]; 
  1054. y[j]=w[j]+h*d[j]/2.0;
  1055.             w[n+j]=y[j];
  1056.         }
  1057.         s=t+h/2.0;
  1058.         FunctionValueTI(s,y,d);
  1059.         for(j=0; j<n; j++)
  1060.         { 
  1061. q[n+j]=d[j];
  1062.             y[j]=w[j]+h*d[j]/2.0;
  1063.             w[n+n+j]=y[j];
  1064.         }
  1065.         FunctionValueTI(s,y,d);
  1066.         for(j=0; j<n; j++) q[n+n+j]=d[j];
  1068.         for(j=0; j<n; j++)
  1069.         { 
  1070. aa=q[n+n+j]-q[n+j];
  1071.             bb=w[n+n+j]-w[n+j];
  1072.             if(-aa*bb*h>0.0)
  1073.             {
  1074. p[j]=-aa/bb; dd=-p[j]*h;
  1075.                 r[j]=exp(dd);
  1076.                 r[n+j]=(r[j]-1.0)/dd;
  1077.                 r[n+n+j]=(r[n+j]-1.0)/dd;
  1078.                 r[3*n+j]=(r[n+n+j]-1.0)/dd;
  1079.             }
  1080.             else p[j]=0.0;
  1081.             if(p[j]<0.0||FloatEqual(p[j],0.0)) g=q[n+n+j];
  1082.             else
  1083.             { 
  1084. g=2.0*(q[n+n+j]-q[j])*r[n+n+j];
  1085.                 g=g+(q[j]-q[n+j])*r[n+j]+q[n+j];
  1086.             }
  1087.             w[3*n+j]=w[j]+g*h;
  1088.             y[j]=w[3*n+j];
  1089.         }
  1090.         s=t+h;
  1091.         FunctionValueTI(s,y,d);
  1092.         for(j=0; j<n; j++) q[3*n+j]=d[j];
  1093.         for(j=0; j<n; j++)
  1094.         { 
  1095. if(p[j]<0.0||FloatEqual(p[j],0.0))
  1096.             { 
  1097. dy=q[j]+2.0*(q[n+j]+q[n+n+j]);
  1098.                 dy=(dy+q[n+n+n+j])*h/6.0;
  1099.             }
  1100.             else
  1101.             { 
  1102. dy=-3.0*(q[j]+p[j]*w[j])+2.0*(q[n+j]+p[j]*w[n+j]);
  1103.                 dy=dy+2.0*(q[n+n+j]+p[j]*w[n+n+j]);
  1104.                 dy=dy-(q[n+n+n+j]+p[j]*w[n+n+n+j]);
  1105.                 dy=dy*r[n+n+j]+q[j]*r[n+j];
  1106.                 dy1=q[j]-q[n+j]-q[n+n+j]+q[n+n+n+j];
  1107.                 dy1=dy1+(w[j]-w[n+j]-w[n+n+j]+w[n+n+n+j])*p[j];
  1108.                 dy=(dy+4.0*dy1*r[n+n+n+j])*h;
  1109.             }
  1110.             y[j]=w[j]+dy;
  1111.             z(j,i)=y[j];
  1112.         }
  1113.     }
  1114. }
  1116. //积分刚性方程组吉尔法
  1117. template <class _Ty>
  1118. int ODE_Gear(_Ty a, _Ty b, _Ty hmin, _Ty hmax, _Ty h, _Ty eps, 
  1119. valarray<_Ty>& y0, int k, valarray<_Ty>& t, matrix<_Ty>& z)
  1120. {
  1121.     int kf(1),jt(0),nn(0),nq(1),m(2),irt1,irt(1),j,nqd,idb;
  1122.     int iw,j1,j2,nt,nqw,l;
  1123.     _Ty aa[7],hw(h),hd,rm,t0(a),td,r,dd,pr1,pr2,pr3,rr;
  1124.     _Ty enq1,enq2,enq3,eup,e,edwn,bnd,r1;
  1126.     _Ty pp[7][3] = 
  1127. {
  1128. {2.0,3.0,1.0}, {4.5,6.0,1.0},
  1129.         {7.333,9.167,0.5}, {10.42,12.5,0.1667},
  1130.         {13.7,15.98,0.04133}, {17.15,1.0,0.008267},
  1131.         {1.0,1.0,1.0}
  1132. };
  1134. int n = y0.size(); //微分方程组中方程的个数，也是未知函数的个数
  1135. valarray<_Ty> d(n), s(10*n), s02(n), ym(n), er(n);
  1136. valarray<_Ty> yy(n), y(8*n), iis(n), jjs(n);
  1137. matrix<_Ty> p(n,n);
  1139.     aa[1]=-1.0;
  1140.     for(int i=0; i<8*n; i++) y[i]=0.0;
  1141.     for(i=0; i<n; i++) 
  1142.     {
  1143. y[i*8]=y0[i];
  1144. yy[i]=y[i*8];
  1145. }
  1146.     FunctionValueG(t0,yy,d);
  1147.     for(i=0; i<n; i++) y[i*8+1]=h*d[i];
  1148.     for(i=0; i<n; i++) ym[i]=1.0;
  1149.   l20:
  1150.     irt=1;
  1151. kf=1; 
  1152. nn=nn+1;
  1153.     t[nn-1]=t0;
  1154.     for(i=0; i<n; i++) z(i,nn-1)=y[i*8];
  1155.     if((t0>b||FloatEqual(t0,b))||(nn==k)) return(kf);
  1156.     
  1157. for(i=0; i<n; i++)
  1158.       for(j=0; j<m; j++) s[i*10+j]=y[i*8+j];
  1159.     hd=hw;
  1160.     if(h!=hd)
  1161.     {
  1162. rm=h/hd;
  1163. irt1=0;
  1164.         rr=Abs(hmin/hd);
  1165.         if(rm<rr) rm=rr;
  1166.         rr=Abs(hmax/hd);
  1167.         if(rm>rr) rm=rr;
  1168. r=1.0; 
  1169. irt1=irt1+1;
  1170.         for(j=1; j<m; j++)
  1171.         { 
  1172. r=r*rm;
  1173.             for(i=0; i<n; i++)
  1174.               y[i*8+j]=s[i*10+j]*r;
  1175.         }
  1176.         h=hd*rm;
  1177.         for(i=0; i<n; i++)
  1178. y[i*8]=s[i*10];
  1179.         idb=m;
  1180.     }
  1181.     nqd=nq;
  1182. td=t0;
  1183. rm=1.0;
  1184.     if(jt>0) goto l80;
  1185.   l60:
  1186.     switch(nq)
  1187.     { 
  1188. case 1: aa[0]=-1.0;
  1189. break;
  1190.         case 2: aa[0]=-2.0/3.0;
  1191. aa[2]=-1.0/3.0; 
  1192. break;
  1193.         case 3: aa[0]=-6.0/11.0; 
  1194. aa[2]=aa[0];
  1195.                 aa[3]=-1.0/11.0; 
  1196. break;
  1197.         case 4: aa[0]=-0.48; 
  1198. aa[2]=-0.7;
  1199. aa[3]=-0.2;
  1200.                 aa[4]=-0.02; 
  1201. break;
  1202.         case 5: aa[0]=-120.0/274.0;
  1203. aa[2]=-225.0/274.0;
  1204.                 aa[3]=-85.0/274.0; 
  1205. aa[4]=-15.0/274.0;
  1206.                 aa[5]=-1.0/274.0; 
  1207. break;
  1208.         case 6: aa[0]=-720.0/1764.0; 
  1209. aa[2]=-1624.0/1764.0;
  1210.                 aa[3]=-735.0/1764.0;
  1211. aa[4]=-175.0/1764.0;
  1212.                 aa[5]=-21.0/1764.0; 
  1213. aa[6]=-1.0/1764.0;
  1214.                 break;
  1215.         default: return(-2);
  1216.     }
  1217.     m=nq+1; 
  1218. idb=m;
  1219.     enq2=0.5/(nq+1.0);
  1220. enq3=0.5/(nq+2.0);
  1221.     enq1=0.5/(nq+0.0);
  1222.     eup=pp[nq-1][1]*eps;
  1223. eup=eup*eup;
  1224.     e=pp[nq-1][0]*eps;
  1225. e=e*e;
  1226.     edwn=pp[nq-1][2]*eps;
  1227. edwn=edwn*edwn;
  1228.     if(edwn==0.0)
  1229.     {
  1230. for(i=0; i<n; i++)
  1231.           for(j=0; j<m; j++)
  1232.             y[i*8+j]=s[i*10+j];
  1233.         h=hd;
  1234. nq=nqd; 
  1235. jt=nq;
  1236. return(-4);
  1237.     }
  1238.     bnd=eps*enq3/(n+0.0);
  1239.     iw=1;
  1240.     if(irt==2)
  1241.     { 
  1242. r1=1.0;
  1243.         for(j=1; j<m; j++)
  1244.         {
  1245. r1=r1*r;
  1246.             for(i=0; i<n; i++)
  1247. y[i*8+j]=y[i*8+j]*r1;
  1248.         }
  1249.         idb=m;
  1250.         for(i=0; i<n; i++)
  1251.           if(ym[i]<Abs(y[i*8]))
  1252. ym[i]=Abs(y[i*8]);
  1253.         jt=nq;
  1254.         goto l20;
  1255.     }
  1256.   l80:
  1257.     t0=t0+h;
  1258.     for(j=2; j<=m; j++)
  1259.       for(j1=j; j1<=m; j1++)
  1260.       { 
  1261.   j2=m-j1+j-1;
  1262.           for(i=0; i<n; i++)
  1263. y[i*8+j2-1]=y[i*8+j2-1]+y[i*8+j2];
  1264.       }
  1265.     for(i=0; i<n; i++) er[i]=0.0;
  1266.     j1=1; 
  1267. nt=1;
  1268.     for(l=0; l<=2; l++)
  1269.     {
  1270. if((j1!=0)&&(nt!=0))
  1271.         {
  1272. for(i=0; i<n; i++) yy[i]=y[i*8];
  1273.             FunctionValueG(t0,yy,d);
  1274.             if(iw>=1)
  1275.             {
  1276. for(i=0; i<n; i++) yy[i]=y[i*8];
  1277.                 JacobiMatrix(t0,yy,p); //计算雅可比矩阵的函数
  1278.                 r=aa[0]*h;
  1279.                 for(i=0; i<n; i++)
  1280.                   for(j=0; j<n; j++)
  1281.                     p(i,j)=p(i,j)*r;
  1282.                 for(i=0; i<n; i++)
  1283.                   p(i,i)=1.0+p(i,i);
  1284.                 iw=-1;
  1285.                 jjs[0]=MatrixInversionGS(p); //全选主元高斯-约当法求矩阵逆
  1286.                 j1=jjs[0];
  1287.             }
  1288.             if(jjs[0]!=0)
  1289.             {
  1290. for(i=0; i<n; i++)
  1291. s02[i]=y[i*8+1]-d[i]*h;
  1292.                 for(i=0; i<n; i++)
  1293.                 { 
  1294. dd=0.0;
  1295.                     for(j=0; j<n; j++)
  1296. dd=dd+s02[j]*p(i,j);
  1297.                     s[i*10+8]=dd;
  1298.                 }
  1299.                 nt=n;
  1300.                 for(i=0; i<n; i++)
  1301.                 { 
  1302. y[i*8]=y[i*8]+aa[0]*s[i*10+8];
  1303.                     y[i*8+1]=y[i*8+1]-s[i*10+8];
  1304.                     er[i]=er[i]+s[i*10+8];
  1305.                     if(Abs(s[i*10+8])<=(bnd*ym[i]))
  1306. nt=nt-1;
  1307.                 }
  1308.             }
  1309.         }
  1310.     }
  1311.     if(nt>0)
  1312.     {
  1313. t0=td;
  1314.         if((h>(hmin*1.00001))||(iw>=0))
  1315.         {
  1316. if(iw!=0) rm=0.25*rm;
  1317.             iw=1;
  1318. irt1=2;
  1319.             rr=Abs(hmin/hd);
  1320.             if(rm<rr) rm=rr;
  1321.             rr=Abs(hmax/hd);
  1322.             if(rm>rr) rm=rr;
  1323.             r=1.0;
  1324.             for(j=1; j<m; j++)
  1325.             {
  1326. r=r*rm;
  1327.                 for(i=0; i<n; i++)
  1328. y[i*8+j]=s[i*10+j]*r;
  1329.             }
  1330.             h=hd*rm;
  1331.             for(i=0; i<n; i++)
  1332.               y[i*8]=s[i*10];
  1333.             idb=m;
  1334.             goto l80;
  1335.         }
  1336.         for(i=0; i<n; i++)
  1337.           for(j=0; j<m; j++)
  1338.             y[i*8+j]=s[i*10+j];
  1339.         h=hd; 
  1340. nq=nqd;
  1341. jt=nq;
  1342. return(-3);
  1343.     }
  1344.     dd=0.0;
  1345.     for(i=0; i<n; i++)
  1346. dd=dd+(er[i]/ym[i])*(er[i]/ym[i]);
  1347.     iw=0;
  1348.     if(dd<=e)
  1349.     { 
  1350. if(m>=3)
  1351.           for(j=2; j<m; j++)
  1352.             for(i=0; i<n; i++)
  1353.               y[i*8+j]=y[i*8+j]+aa[j]*er[i];
  1354.         kf=1;
  1355. hw=h;
  1356.         if(idb>1)
  1357.         {
  1358. idb=idb-1;
  1359.             if(idb<=1)
  1360.               for(i=0; i<n; i++)
  1361.                 s[i*10+9]=er[i];
  1362.             for(i=0; i<n; i++)
  1363.       if(ym[i]<Abs(y[i*8])) ym[i]=Abs(y[i*8]);
  1364.             jt=nq;
  1365.             goto l20;
  1366.         }
  1367.     }
  1368.     if(dd>e)
  1369.     { 
  1370. kf=kf-2;
  1371.         if(h<=(hmin*1.00001))
  1372.         {
  1373.             hw=h;
  1374. jt=nq; 
  1375. return(-1);
  1376.         }
  1377.         t0=td;
  1378.         if(kf<=-5)
  1379.         { 
  1380. if(nq==1)
  1381.             {
  1382. for(i=0; i<n; i++)
  1383.                   for(j=0; j<m; j++)
  1384.                     y[i*8+j]=s[i*10+j];
  1385.                 h=hd;
  1386. nq=nqd;
  1387. jt=nq;
  1388. return(-4);
  1389.             }
  1390.             for(i=0; i<n; i++) yy[i]=y[i*8];
  1391.             FunctionValueG(t0,yy,d);
  1392.             r=h/hd;
  1393.             for(i=0; i<n; i++)
  1394.             { 
  1395. y[i*8]=s[i*10];
  1396.                 s[i*10+1]=hd*d[i];
  1397.                 y[i*8+1]=s[i*10+1]*r;
  1398.             }
  1399.             nq=1;
  1400. kf=1;
  1401. goto l60;
  1402.         }
  1403.     }
  1404.     pr2=log(dd/e);
  1405. pr2=enq2*pr2;
  1406. pr2=exp(pr2);
  1407.     pr2=1.2*pr2;
  1408.     pr3=1.0e+20;
  1409.     if(nq<7)
  1410.       if(kf>-1)
  1411.       {
  1412.   dd=0.0;
  1413.           for(i=0; i<n; i++)
  1414.           {
  1415.   pr3=(er[i]-s[i*10+9])/ym[i];
  1416.               dd=dd+pr3*pr3;
  1417.           }
  1418.           pr3=log(dd/eup);
  1419.   pr3=enq3*pr3;
  1420.           pr3=exp(pr3);
  1421.   pr3=1.4*pr3;
  1422.       }
  1423.     pr1=1.0e+20;
  1424.     if(nq>1)
  1425.     { 
  1426. dd=0.0;
  1427.         for(i=0; i<n; i++)
  1428.         {
  1429. pr1=y[i*8+m-1]/ym[i];
  1430.             dd=dd+pr1*pr1;
  1431.         }
  1432.         pr1=log(dd/edwn); 
  1433. pr1=enq1*pr1;
  1434.         pr1=exp(pr1);
  1435. pr1=1.3*pr1;
  1436.     }
  1437.     if(pr2<=pr3)
  1438.     { 
  1439. if(pr2>pr1)
  1440. {
  1441. r=1.0e+04;
  1442.             if(pr1>1.0e-04) r=1.0/pr1;
  1443.             nqw=nq-1;
  1444.         }
  1445.         else
  1446.         { 
  1447. nqw=nq;
  1448. r=1.0e+04;
  1449.             if(pr2>1.0e-04) r=1.0/pr2;
  1450.         }
  1451.     }
  1452.     else
  1453.     { 
  1454. if(pr3<pr1)
  1455.         { 
  1456. r=1.0e+04;
  1457.             if(pr3>1.0e-04) r=1.0/pr3;
  1458.             nqw=nq+1;
  1459.         }
  1460.         else
  1461.         { 
  1462. r=1.0e+04;
  1463.             if(pr1>1.0e-04) r=1.0/pr1;
  1464.             nqw=nq-1;
  1465.         }
  1466.     }
  1467.     idb=10;
  1468.     if(kf==1)
  1469.       if(r<1.1)
  1470.       { 
  1471.   for(i=0; i<n; i++)
  1472.             if(ym[i]<Abs(y[i*8])) ym[i]=Abs(y[i*8]);
  1473.           jt=nq; goto l20;
  1474.       }
  1475.     if(nqw>nq)
  1476.       for(i=0; i<n; i++)
  1477.         y[i*8+nqw]=er[i]*aa[m-1]/(m+0.0);
  1478.     m=nqw+1;
  1479.     if(kf==1)
  1480.     { 
  1481. irt=2; 
  1482. rr=hmax/Abs(h);
  1483.         if(r>rr) r=rr;
  1484.         h=h*r; 
  1485. hw=h;
  1486.         if(nq==nqw)
  1487.         {
  1488. r1=1.0;
  1489.             for(j=1; j<m; j++)
  1490.             {
  1491. r1=r1*r;
  1492.                 for(i=0; i<n; i++)
  1493.                   y[i*8+j]=y[i*8+j]*r1;
  1494.             }
  1495.             idb=m;
  1496.             for(i=0; i<n; i++)
  1497.               if(ym[i]<Abs(y[i*8])) ym[i]=Abs(y[i*8]);
  1498.             jt=nq;
  1499. goto l20;
  1500.         }
  1501.         nq=nqw;
  1502.         goto l60;
  1503.     }
  1504.     rm=rm*r;
  1505. irt1=3;
  1506.     rr=Abs(hmin/hd);
  1507.     if(rm<rr) rm=rr;
  1508.     rr=Abs(hmax/hd);
  1509.     if(rm>rr) rm=rr;
  1510.     r=1.0;
  1511.     for(j=1; j<m; j++)
  1512.     {
  1513. r=r*rm;
  1514.         for(i=0; i<n; i++)
  1515.           y[i*8+j]=s[i*10+j]*r;
  1516.     }
  1517.     h=hd*rm;
  1518.     for(i=0; i<n; i++)
  1519.       y[i*8]=s[i*10];
  1520.     idb=m;
  1521.     if(nqw==nq) goto l80;
  1522.     nq=nqw; goto l60;
  1523. }
  1525. //二阶微分方程边值问题数值解法
  1526. template <class _Ty>
  1527. void ODE_LinearBoundaryValude(_Ty a, _Ty b, _Ty ya, _Ty yb, int n, 
  1528. valarray<_Ty>& y)
  1529. {
  1530.     int i,k,nn,m1;
  1531.     _Ty h,x;
  1532. valarray<_Ty> g(6*n), d(2*n), z(4);
  1534.     h=(b-a)/(n-1.0);
  1535. nn=2*n-1;
  1536.     g[0]=1.0;
  1537. g[1]=0.0;
  1538.     y[0]=ya; 
  1539. y[n-1]=yb;
  1540.     g[3*n-3]=1.0;
  1541. g[3*n-4]=0.0;
  1542.     for(i=2; i<n; i++)
  1543.     { 
  1544. x=a+(i-1)*h;
  1545.         FunctionValueLBV(x,z);
  1546.         k=3*(i-1)-1;
  1547.         g[k]=z[0]-h*z[1]/2.0;
  1548.         g[k+1]=h*h*z[2]-2.0*z[0];
  1549.         g[k+2]=z[0]+h*z[1]/2.0;
  1550.         y[i-1]=h*h*z[3];
  1551.     }
  1552.     m1=3*n-2;
  1554. valarray<_Ty> gg(m1);
  1555. for(i=0;i<m1;i++) gg[i] = g[i];
  1556. LE_TridiagonalEquationGauss(gg,y); //求解三对角方程组
  1558. for(i=0;i<m1;i++) g[i] = gg[i];
  1560.     h=h/2.0;
  1561.     g[0]=1.0; 
  1562. g[1]=0.0;
  1563.     d[0]=ya; 
  1564. d[nn-1]=yb;
  1565.     g[3*nn-3]=1.0;
  1566. g[3*nn-4]=0.0;
  1567.     for(i=2; i<nn; i++)
  1568.     { 
  1569. x=a+(i-1)*h;
  1570.         FunctionValueLBV(x,z);
  1571.         k=3*(i-1)-1;
  1572.         g[k]=z[0]-h*z[1]/2.0;
  1573.         g[k+1]=h*h*z[2]-2.0*z[0];
  1574.         g[k+2]=z[0]+h*z[1]/2.0;
  1575.         d[i-1]=h*h*z[3];
  1576.     }
  1577.     m1=3*nn-2;
  1579. valarray<_Ty> ggg(m1), ddd(nn);
  1580. for(i=0;i<m1;i++) ggg[i] = g[i];
  1581. for(i=0;i<nn;i++) ddd[i] = d[i];
  1582. LE_TridiagonalEquationGauss(ggg,ddd); //求解三对角方程组
  1584.     for(i=0;i<m1;i++) g[i] = ggg[i];
  1585. for(i=0;i<nn;i++) d[i] = ddd[i];
  1587. for(i=2; i<n; i++)
  1588.     { 
  1589. k=i+i-1;
  1590.         y[i-1]=(4.0*d[k-1]-y[i-1])/3.0;
  1591.     }
  1592. }
  1594. #endif // _ORDINARYDIFFERENTIALEGUATION_INL