1  /*************************************************************************


2  Copyright (c) Sergey Bochkanov (ALGLIB project).


3 


4  >>> SOURCE LICENSE >>>


5  This program is free software; you can redistribute it and/or modify


6  it under the terms of the GNU General Public License as published by


7  the Free Software Foundation (www.fsf.org); either version 2 of the


8  License, or (at your option) any later version.


9 


10  This program is distributed in the hope that it will be useful,


11  but WITHOUT ANY WARRANTY; without even the implied warranty of


12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


13  GNU General Public License for more details.


14 


15  A copy of the GNU General Public License is available at


16  http://www.fsf.org/licensing/licenses


17  >>> END OF LICENSE >>>


18  *************************************************************************/


19  #pragma warning disable 162


20  #pragma warning disable 219


21  using System;


22 


23  public partial class alglib


24  {


25 


26 


27  /*************************************************************************


28 


29  *************************************************************************/


30  public class odesolverstate


31  {


32  //


33  // Public declarations


34  //


35  public bool needdy { get { return _innerobj.needdy; } set { _innerobj.needdy = value; } }


36  public double[] y { get { return _innerobj.y; } }


37  public double[] dy { get { return _innerobj.dy; } }


38  public double x { get { return _innerobj.x; } set { _innerobj.x = value; } }


39 


40  public odesolverstate()


41  {


42  _innerobj = new odesolver.odesolverstate();


43  }


44 


45  //


46  // Although some of declarations below are public, you should not use them


47  // They are intended for internal use only


48  //


49  private odesolver.odesolverstate _innerobj;


50  public odesolver.odesolverstate innerobj { get { return _innerobj; } }


51  public odesolverstate(odesolver.odesolverstate obj)


52  {


53  _innerobj = obj;


54  }


55  }


56 


57 


58  /*************************************************************************


59 


60  *************************************************************************/


61  public class odesolverreport


62  {


63  //


64  // Public declarations


65  //


66  public int nfev { get { return _innerobj.nfev; } set { _innerobj.nfev = value; } }


67  public int terminationtype { get { return _innerobj.terminationtype; } set { _innerobj.terminationtype = value; } }


68 


69  public odesolverreport()


70  {


71  _innerobj = new odesolver.odesolverreport();


72  }


73 


74  //


75  // Although some of declarations below are public, you should not use them


76  // They are intended for internal use only


77  //


78  private odesolver.odesolverreport _innerobj;


79  public odesolver.odesolverreport innerobj { get { return _innerobj; } }


80  public odesolverreport(odesolver.odesolverreport obj)


81  {


82  _innerobj = obj;


83  }


84  }


85 


86  /*************************************************************************


87  CashKarp adaptive ODE solver.


88 


89  This subroutine solves ODE Y'=f(Y,x) with initial conditions Y(xs)=Ys


90  (here Y may be single variable or vector of N variables).


91 


92  INPUT PARAMETERS:


93  Y  initial conditions, array[0..N1].


94  contains values of Y[] at X[0]


95  N  system size


96  X  points at which Y should be tabulated, array[0..M1]


97  integrations starts at X[0], ends at X[M1], intermediate


98  values at X[i] are returned too.


99  SHOULD BE ORDERED BY ASCENDING OR BY DESCENDING!!!!


100  M  number of intermediate points + first point + last point:


101  * M>2 means that you need both Y(X[M1]) and M2 values at


102  intermediate points


103  * M=2 means that you want just to integrate from X[0] to


104  X[1] and don't interested in intermediate values.


105  * M=1 means that you don't want to integrate :)


106  it is degenerate case, but it will be handled correctly.


107  * M<1 means error


108  Eps  tolerance (absolute/relative error on each step will be


109  less than Eps). When passing:


110  * Eps>0, it means desired ABSOLUTE error


111  * Eps<0, it means desired RELATIVE error. Relative errors


112  are calculated with respect to maximum values of Y seen


113  so far. Be careful to use this criterion when starting


114  from Y[] that are close to zero.


115  H  initial step lenth, it will be adjusted automatically


116  after the first step. If H=0, step will be selected


117  automatically (usualy it will be equal to 0.001 of


118  min(x[i]x[j])).


119 


120  OUTPUT PARAMETERS


121  State  structure which stores algorithm state between subsequent


122  calls of OdeSolverIteration. Used for reverse communication.


123  This structure should be passed to the OdeSolverIteration


124  subroutine.


125 


126  SEE ALSO


127  AutoGKSmoothW, AutoGKSingular, AutoGKIteration, AutoGKResults.


128 


129 


130   ALGLIB 


131  Copyright 01.09.2009 by Bochkanov Sergey


132  *************************************************************************/


133  public static void odesolverrkck(double[] y, int n, double[] x, int m, double eps, double h, out odesolverstate state)


134  {


135  state = new odesolverstate();


136  odesolver.odesolverrkck(y, n, x, m, eps, h, state.innerobj);


137  return;


138  }


139  public static void odesolverrkck(double[] y, double[] x, double eps, double h, out odesolverstate state)


140  {


141  int n;


142  int m;


143 


144  state = new odesolverstate();


145  n = ap.len(y);


146  m = ap.len(x);


147  odesolver.odesolverrkck(y, n, x, m, eps, h, state.innerobj);


148 


149  return;


150  }


151 


152  /*************************************************************************


153  This function provides reverse communication interface


154  Reverse communication interface is not documented or recommended to use.


155  See below for functions which provide better documented API


156  *************************************************************************/


157  public static bool odesolveriteration(odesolverstate state)


158  {


159 


160  bool result = odesolver.odesolveriteration(state.innerobj);


161  return result;


162  }


163  /*************************************************************************


164  This function is used to launcn iterations of ODE solver


165 


166  It accepts following parameters:


167  diff  callback which calculates dy/dx for given y and x


168  obj  optional object which is passed to diff; can be NULL


169 


170 


171   ALGLIB 


172  Copyright 01.09.2009 by Bochkanov Sergey


173 


174  *************************************************************************/


175  public static void odesolversolve(odesolverstate state, ndimensional_ode_rp diff, object obj)


176  {


177  if( diff==null )


178  throw new alglibexception("ALGLIB: error in 'odesolversolve()' (diff is null)");


179  while( alglib.odesolveriteration(state) )


180  {


181  if( state.needdy )


182  {


183  diff(state.innerobj.y, state.innerobj.x, state.innerobj.dy, obj);


184  continue;


185  }


186  throw new alglibexception("ALGLIB: unexpected error in 'odesolversolve'");


187  }


188  }


189 


190 


191 


192  /*************************************************************************


193  ODE solver results


194 


195  Called after OdeSolverIteration returned False.


196 


197  INPUT PARAMETERS:


198  State  algorithm state (used by OdeSolverIteration).


199 


200  OUTPUT PARAMETERS:


201  M  number of tabulated values, M>=1


202  XTbl  array[0..M1], values of X


203  YTbl  array[0..M1,0..N1], values of Y in X[i]


204  Rep  solver report:


205  * Rep.TerminationType completetion code:


206  * 2 X is not ordered by ascending/descending or


207  there are nondistinct X[], i.e. X[i]=X[i+1]


208  * 1 incorrect parameters were specified


209  * 1 task has been solved


210  * Rep.NFEV contains number of function calculations


211 


212   ALGLIB 


213  Copyright 01.09.2009 by Bochkanov Sergey


214  *************************************************************************/


215  public static void odesolverresults(odesolverstate state, out int m, out double[] xtbl, out double[,] ytbl, out odesolverreport rep)


216  {


217  m = 0;


218  xtbl = new double[0];


219  ytbl = new double[0,0];


220  rep = new odesolverreport();


221  odesolver.odesolverresults(state.innerobj, ref m, ref xtbl, ref ytbl, rep.innerobj);


222  return;


223  }


224 


225  }


226  public partial class alglib


227  {


228  public class odesolver


229  {


230  public class odesolverstate


231  {


232  public int n;


233  public int m;


234  public double xscale;


235  public double h;


236  public double eps;


237  public bool fraceps;


238  public double[] yc;


239  public double[] escale;


240  public double[] xg;


241  public int solvertype;


242  public bool needdy;


243  public double x;


244  public double[] y;


245  public double[] dy;


246  public double[,] ytbl;


247  public int repterminationtype;


248  public int repnfev;


249  public double[] yn;


250  public double[] yns;


251  public double[] rka;


252  public double[] rkc;


253  public double[] rkcs;


254  public double[,] rkb;


255  public double[,] rkk;


256  public rcommstate rstate;


257  public odesolverstate()


258  {


259  yc = new double[0];


260  escale = new double[0];


261  xg = new double[0];


262  y = new double[0];


263  dy = new double[0];


264  ytbl = new double[0,0];


265  yn = new double[0];


266  yns = new double[0];


267  rka = new double[0];


268  rkc = new double[0];


269  rkcs = new double[0];


270  rkb = new double[0,0];


271  rkk = new double[0,0];


272  rstate = new rcommstate();


273  }


274  };


275 


276 


277  public class odesolverreport


278  {


279  public int nfev;


280  public int terminationtype;


281  };


282 


283 


284 


285 


286  public const double odesolvermaxgrow = 3.0;


287  public const double odesolvermaxshrink = 10.0;


288 


289 


290  /*************************************************************************


291  CashKarp adaptive ODE solver.


292 


293  This subroutine solves ODE Y'=f(Y,x) with initial conditions Y(xs)=Ys


294  (here Y may be single variable or vector of N variables).


295 


296  INPUT PARAMETERS:


297  Y  initial conditions, array[0..N1].


298  contains values of Y[] at X[0]


299  N  system size


300  X  points at which Y should be tabulated, array[0..M1]


301  integrations starts at X[0], ends at X[M1], intermediate


302  values at X[i] are returned too.


303  SHOULD BE ORDERED BY ASCENDING OR BY DESCENDING!!!!


304  M  number of intermediate points + first point + last point:


305  * M>2 means that you need both Y(X[M1]) and M2 values at


306  intermediate points


307  * M=2 means that you want just to integrate from X[0] to


308  X[1] and don't interested in intermediate values.


309  * M=1 means that you don't want to integrate :)


310  it is degenerate case, but it will be handled correctly.


311  * M<1 means error


312  Eps  tolerance (absolute/relative error on each step will be


313  less than Eps). When passing:


314  * Eps>0, it means desired ABSOLUTE error


315  * Eps<0, it means desired RELATIVE error. Relative errors


316  are calculated with respect to maximum values of Y seen


317  so far. Be careful to use this criterion when starting


318  from Y[] that are close to zero.


319  H  initial step lenth, it will be adjusted automatically


320  after the first step. If H=0, step will be selected


321  automatically (usualy it will be equal to 0.001 of


322  min(x[i]x[j])).


323 


324  OUTPUT PARAMETERS


325  State  structure which stores algorithm state between subsequent


326  calls of OdeSolverIteration. Used for reverse communication.


327  This structure should be passed to the OdeSolverIteration


328  subroutine.


329 


330  SEE ALSO


331  AutoGKSmoothW, AutoGKSingular, AutoGKIteration, AutoGKResults.


332 


333 


334   ALGLIB 


335  Copyright 01.09.2009 by Bochkanov Sergey


336  *************************************************************************/


337  public static void odesolverrkck(double[] y,


338  int n,


339  double[] x,


340  int m,


341  double eps,


342  double h,


343  odesolverstate state)


344  {


345  alglib.ap.assert(n>=1, "ODESolverRKCK: N<1!");


346  alglib.ap.assert(m>=1, "ODESolverRKCK: M<1!");


347  alglib.ap.assert(alglib.ap.len(y)>=n, "ODESolverRKCK: Length(Y)<N!");


348  alglib.ap.assert(alglib.ap.len(x)>=m, "ODESolverRKCK: Length(X)<M!");


349  alglib.ap.assert(apserv.isfinitevector(y, n), "ODESolverRKCK: Y contains infinite or NaN values!");


350  alglib.ap.assert(apserv.isfinitevector(x, m), "ODESolverRKCK: Y contains infinite or NaN values!");


351  alglib.ap.assert(math.isfinite(eps), "ODESolverRKCK: Eps is not finite!");


352  alglib.ap.assert((double)(eps)!=(double)(0), "ODESolverRKCK: Eps is zero!");


353  alglib.ap.assert(math.isfinite(h), "ODESolverRKCK: H is not finite!");


354  odesolverinit(0, y, n, x, m, eps, h, state);


355  }


356 


357 


358  /*************************************************************************


359 


360   ALGLIB 


361  Copyright 01.09.2009 by Bochkanov Sergey


362  *************************************************************************/


363  public static bool odesolveriteration(odesolverstate state)


364  {


365  bool result = new bool();


366  int n = 0;


367  int m = 0;


368  int i = 0;


369  int j = 0;


370  int k = 0;


371  double xc = 0;


372  double v = 0;


373  double h = 0;


374  double h2 = 0;


375  bool gridpoint = new bool();


376  double err = 0;


377  double maxgrowpow = 0;


378  int klimit = 0;


379  int i_ = 0;


380 


381 


382  //


383  // Reverse communication preparations


384  // I know it looks ugly, but it works the same way


385  // anywhere from C++ to Python.


386  //


387  // This code initializes locals by:


388  // * random values determined during code


389  // generation  on first subroutine call


390  // * values from previous call  on subsequent calls


391  //


392  if( state.rstate.stage>=0 )


393  {


394  n = state.rstate.ia[0];


395  m = state.rstate.ia[1];


396  i = state.rstate.ia[2];


397  j = state.rstate.ia[3];


398  k = state.rstate.ia[4];


399  klimit = state.rstate.ia[5];


400  gridpoint = state.rstate.ba[0];


401  xc = state.rstate.ra[0];


402  v = state.rstate.ra[1];


403  h = state.rstate.ra[2];


404  h2 = state.rstate.ra[3];


405  err = state.rstate.ra[4];


406  maxgrowpow = state.rstate.ra[5];


407  }


408  else


409  {


410  n = 983;


411  m = 989;


412  i = 834;


413  j = 900;


414  k = 287;


415  klimit = 364;


416  gridpoint = false;


417  xc = 338;


418  v = 686;


419  h = 912;


420  h2 = 585;


421  err = 497;


422  maxgrowpow = 271;


423  }


424  if( state.rstate.stage==0 )


425  {


426  goto lbl_0;


427  }


428 


429  //


430  // Routine body


431  //


432 


433  //


434  // prepare


435  //


436  if( state.repterminationtype!=0 )


437  {


438  result = false;


439  return result;


440  }


441  n = state.n;


442  m = state.m;


443  h = state.h;


444  maxgrowpow = Math.Pow(odesolvermaxgrow, 5);


445  state.repnfev = 0;


446 


447  //


448  // some preliminary checks for internal errors


449  // after this we assume that H>0 and M>1


450  //


451  alglib.ap.assert((double)(state.h)>(double)(0), "ODESolver: internal error");


452  alglib.ap.assert(m>1, "ODESolverIteration: internal error");


453 


454  //


455  // choose solver


456  //


457  if( state.solvertype!=0 )


458  {


459  goto lbl_1;


460  }


461 


462  //


463  // CaskKarp solver


464  // Prepare coefficients table.


465  // Check it for errors


466  //


467  state.rka = new double[6];


468  state.rka[0] = 0;


469  state.rka[1] = (double)1/(double)5;


470  state.rka[2] = (double)3/(double)10;


471  state.rka[3] = (double)3/(double)5;


472  state.rka[4] = 1;


473  state.rka[5] = (double)7/(double)8;


474  state.rkb = new double[6, 5];


475  state.rkb[1,0] = (double)1/(double)5;


476  state.rkb[2,0] = (double)3/(double)40;


477  state.rkb[2,1] = (double)9/(double)40;


478  state.rkb[3,0] = (double)3/(double)10;


479  state.rkb[3,1] = ((double)9/(double)10);


480  state.rkb[3,2] = (double)6/(double)5;


481  state.rkb[4,0] = ((double)11/(double)54);


482  state.rkb[4,1] = (double)5/(double)2;


483  state.rkb[4,2] = ((double)70/(double)27);


484  state.rkb[4,3] = (double)35/(double)27;


485  state.rkb[5,0] = (double)1631/(double)55296;


486  state.rkb[5,1] = (double)175/(double)512;


487  state.rkb[5,2] = (double)575/(double)13824;


488  state.rkb[5,3] = (double)44275/(double)110592;


489  state.rkb[5,4] = (double)253/(double)4096;


490  state.rkc = new double[6];


491  state.rkc[0] = (double)37/(double)378;


492  state.rkc[1] = 0;


493  state.rkc[2] = (double)250/(double)621;


494  state.rkc[3] = (double)125/(double)594;


495  state.rkc[4] = 0;


496  state.rkc[5] = (double)512/(double)1771;


497  state.rkcs = new double[6];


498  state.rkcs[0] = (double)2825/(double)27648;


499  state.rkcs[1] = 0;


500  state.rkcs[2] = (double)18575/(double)48384;


501  state.rkcs[3] = (double)13525/(double)55296;


502  state.rkcs[4] = (double)277/(double)14336;


503  state.rkcs[5] = (double)1/(double)4;


504  state.rkk = new double[6, n];


505 


506  //


507  // Main cycle consists of two iterations:


508  // * outer where we travel from X[i1] to X[i]


509  // * inner where we travel inside [X[i1],X[i]]


510  //


511  state.ytbl = new double[m, n];


512  state.escale = new double[n];


513  state.yn = new double[n];


514  state.yns = new double[n];


515  xc = state.xg[0];


516  for(i_=0; i_<=n1;i_++)


517  {


518  state.ytbl[0,i_] = state.yc[i_];


519  }


520  for(j=0; j<=n1; j++)


521  {


522  state.escale[j] = 0;


523  }


524  i = 1;


525  lbl_3:


526  if( i>m1 )


527  {


528  goto lbl_5;


529  }


530 


531  //


532  // begin inner iteration


533  //


534  lbl_6:


535  if( false )


536  {


537  goto lbl_7;


538  }


539 


540  //


541  // truncate step if needed (beyond right boundary).


542  // determine should we store X or not


543  //


544  if( (double)(xc+h)>=(double)(state.xg[i]) )


545  {


546  h = state.xg[i]xc;


547  gridpoint = true;


548  }


549  else


550  {


551  gridpoint = false;


552  }


553 


554  //


555  // Update error scale maximums


556  //


557  // These maximums are initialized by zeros,


558  // then updated every iterations.


559  //


560  for(j=0; j<=n1; j++)


561  {


562  state.escale[j] = Math.Max(state.escale[j], Math.Abs(state.yc[j]));


563  }


564 


565  //


566  // make one step:


567  // 1. calculate all info needed to do step


568  // 2. update errors scale maximums using values/derivatives


569  // obtained during (1)


570  //


571  // Take into account that we use scaling of X to reduce task


572  // to the form where x[0] < x[1] < ... < x[n1]. So X is


573  // replaced by x=xscale*t, and dy/dx=f(y,x) is replaced


574  // by dy/dt=xscale*f(y,xscale*t).


575  //


576  for(i_=0; i_<=n1;i_++)


577  {


578  state.yn[i_] = state.yc[i_];


579  }


580  for(i_=0; i_<=n1;i_++)


581  {


582  state.yns[i_] = state.yc[i_];


583  }


584  k = 0;


585  lbl_8:


586  if( k>5 )


587  {


588  goto lbl_10;


589  }


590 


591  //


592  // prepare data for the next update of YN/YNS


593  //


594  state.x = state.xscale*(xc+state.rka[k]*h);


595  for(i_=0; i_<=n1;i_++)


596  {


597  state.y[i_] = state.yc[i_];


598  }


599  for(j=0; j<=k1; j++)


600  {


601  v = state.rkb[k,j];


602  for(i_=0; i_<=n1;i_++)


603  {


604  state.y[i_] = state.y[i_] + v*state.rkk[j,i_];


605  }


606  }


607  state.needdy = true;


608  state.rstate.stage = 0;


609  goto lbl_rcomm;


610  lbl_0:


611  state.needdy = false;


612  state.repnfev = state.repnfev+1;


613  v = h*state.xscale;


614  for(i_=0; i_<=n1;i_++)


615  {


616  state.rkk[k,i_] = v*state.dy[i_];


617  }


618 


619  //


620  // update YN/YNS


621  //


622  v = state.rkc[k];


623  for(i_=0; i_<=n1;i_++)


624  {


625  state.yn[i_] = state.yn[i_] + v*state.rkk[k,i_];


626  }


627  v = state.rkcs[k];


628  for(i_=0; i_<=n1;i_++)


629  {


630  state.yns[i_] = state.yns[i_] + v*state.rkk[k,i_];


631  }


632  k = k+1;


633  goto lbl_8;


634  lbl_10:


635 


636  //


637  // estimate error


638  //


639  err = 0;


640  for(j=0; j<=n1; j++)


641  {


642  if( !state.fraceps )


643  {


644 


645  //


646  // absolute error is estimated


647  //


648  err = Math.Max(err, Math.Abs(state.yn[j]state.yns[j]));


649  }


650  else


651  {


652 


653  //


654  // Relative error is estimated


655  //


656  v = state.escale[j];


657  if( (double)(v)==(double)(0) )


658  {


659  v = 1;


660  }


661  err = Math.Max(err, Math.Abs(state.yn[j]state.yns[j])/v);


662  }


663  }


664 


665  //


666  // calculate new step, restart if necessary


667  //


668  if( (double)(maxgrowpow*err)<=(double)(state.eps) )


669  {


670  h2 = odesolvermaxgrow*h;


671  }


672  else


673  {


674  h2 = h*Math.Pow(state.eps/err, 0.2);


675  }


676  if( (double)(h2)<(double)(h/odesolvermaxshrink) )


677  {


678  h2 = h/odesolvermaxshrink;


679  }


680  if( (double)(err)>(double)(state.eps) )


681  {


682  h = h2;


683  goto lbl_6;


684  }


685 


686  //


687  // advance position


688  //


689  xc = xc+h;


690  for(i_=0; i_<=n1;i_++)


691  {


692  state.yc[i_] = state.yn[i_];


693  }


694 


695  //


696  // update H


697  //


698  h = h2;


699 


700  //


701  // break on grid point


702  //


703  if( gridpoint )


704  {


705  goto lbl_7;


706  }


707  goto lbl_6;


708  lbl_7:


709 


710  //


711  // save result


712  //


713  for(i_=0; i_<=n1;i_++)


714  {


715  state.ytbl[i,i_] = state.yc[i_];


716  }


717  i = i+1;


718  goto lbl_3;


719  lbl_5:


720  state.repterminationtype = 1;


721  result = false;


722  return result;


723  lbl_1:


724  result = false;


725  return result;


726 


727  //


728  // Saving state


729  //


730  lbl_rcomm:


731  result = true;


732  state.rstate.ia[0] = n;


733  state.rstate.ia[1] = m;


734  state.rstate.ia[2] = i;


735  state.rstate.ia[3] = j;


736  state.rstate.ia[4] = k;


737  state.rstate.ia[5] = klimit;


738  state.rstate.ba[0] = gridpoint;


739  state.rstate.ra[0] = xc;


740  state.rstate.ra[1] = v;


741  state.rstate.ra[2] = h;


742  state.rstate.ra[3] = h2;


743  state.rstate.ra[4] = err;


744  state.rstate.ra[5] = maxgrowpow;


745  return result;


746  }


747 


748 


749  /*************************************************************************


750  ODE solver results


751 


752  Called after OdeSolverIteration returned False.


753 


754  INPUT PARAMETERS:


755  State  algorithm state (used by OdeSolverIteration).


756 


757  OUTPUT PARAMETERS:


758  M  number of tabulated values, M>=1


759  XTbl  array[0..M1], values of X


760  YTbl  array[0..M1,0..N1], values of Y in X[i]


761  Rep  solver report:


762  * Rep.TerminationType completetion code:


763  * 2 X is not ordered by ascending/descending or


764  there are nondistinct X[], i.e. X[i]=X[i+1]


765  * 1 incorrect parameters were specified


766  * 1 task has been solved


767  * Rep.NFEV contains number of function calculations


768 


769   ALGLIB 


770  Copyright 01.09.2009 by Bochkanov Sergey


771  *************************************************************************/


772  public static void odesolverresults(odesolverstate state,


773  ref int m,


774  ref double[] xtbl,


775  ref double[,] ytbl,


776  odesolverreport rep)


777  {


778  double v = 0;


779  int i = 0;


780  int i_ = 0;


781 


782  m = 0;


783  xtbl = new double[0];


784  ytbl = new double[0,0];


785 


786  rep.terminationtype = state.repterminationtype;


787  if( rep.terminationtype>0 )


788  {


789  m = state.m;


790  rep.nfev = state.repnfev;


791  xtbl = new double[state.m];


792  v = state.xscale;


793  for(i_=0; i_<=state.m1;i_++)


794  {


795  xtbl[i_] = v*state.xg[i_];


796  }


797  ytbl = new double[state.m, state.n];


798  for(i=0; i<=state.m1; i++)


799  {


800  for(i_=0; i_<=state.n1;i_++)


801  {


802  ytbl[i,i_] = state.ytbl[i,i_];


803  }


804  }


805  }


806  else


807  {


808  rep.nfev = 0;


809  }


810  }


811 


812 


813  /*************************************************************************


814  Internal initialization subroutine


815  *************************************************************************/


816  private static void odesolverinit(int solvertype,


817  double[] y,


818  int n,


819  double[] x,


820  int m,


821  double eps,


822  double h,


823  odesolverstate state)


824  {


825  int i = 0;


826  double v = 0;


827  int i_ = 0;


828 


829 


830  //


831  // Prepare RComm


832  //


833  state.rstate.ia = new int[5+1];


834  state.rstate.ba = new bool[0+1];


835  state.rstate.ra = new double[5+1];


836  state.rstate.stage = 1;


837  state.needdy = false;


838 


839  //


840  // check parameters.


841  //


842  if( (n<=0  m<1)  (double)(eps)==(double)(0) )


843  {


844  state.repterminationtype = 1;


845  return;


846  }


847  if( (double)(h)<(double)(0) )


848  {


849  h = h;


850  }


851 


852  //


853  // quick exit if necessary.


854  // after this block we assume that M>1


855  //


856  if( m==1 )


857  {


858  state.repnfev = 0;


859  state.repterminationtype = 1;


860  state.ytbl = new double[1, n];


861  for(i_=0; i_<=n1;i_++)


862  {


863  state.ytbl[0,i_] = y[i_];


864  }


865  state.xg = new double[m];


866  for(i_=0; i_<=m1;i_++)


867  {


868  state.xg[i_] = x[i_];


869  }


870  return;


871  }


872 


873  //


874  // check again: correct order of X[]


875  //


876  if( (double)(x[1])==(double)(x[0]) )


877  {


878  state.repterminationtype = 2;


879  return;


880  }


881  for(i=1; i<=m1; i++)


882  {


883  if( ((double)(x[1])>(double)(x[0]) && (double)(x[i])<=(double)(x[i1]))  ((double)(x[1])<(double)(x[0]) && (double)(x[i])>=(double)(x[i1])) )


884  {


885  state.repterminationtype = 2;


886  return;


887  }


888  }


889 


890  //


891  // autoselect H if necessary


892  //


893  if( (double)(h)==(double)(0) )


894  {


895  v = Math.Abs(x[1]x[0]);


896  for(i=2; i<=m1; i++)


897  {


898  v = Math.Min(v, Math.Abs(x[i]x[i1]));


899  }


900  h = 0.001*v;


901  }


902 


903  //


904  // store parameters


905  //


906  state.n = n;


907  state.m = m;


908  state.h = h;


909  state.eps = Math.Abs(eps);


910  state.fraceps = (double)(eps)<(double)(0);


911  state.xg = new double[m];


912  for(i_=0; i_<=m1;i_++)


913  {


914  state.xg[i_] = x[i_];


915  }


916  if( (double)(x[1])>(double)(x[0]) )


917  {


918  state.xscale = 1;


919  }


920  else


921  {


922  state.xscale = 1;


923  for(i_=0; i_<=m1;i_++)


924  {


925  state.xg[i_] = 1*state.xg[i_];


926  }


927  }


928  state.yc = new double[n];


929  for(i_=0; i_<=n1;i_++)


930  {


931  state.yc[i_] = y[i_];


932  }


933  state.solvertype = solvertype;


934  state.repterminationtype = 0;


935 


936  //


937  // Allocate arrays


938  //


939  state.y = new double[n];


940  state.dy = new double[n];


941  }


942 


943 


944  }


945  }


946 

