home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1994 #1 / monster.zip / monster / PROG_C / C_LSODE.ZIP / TESTE.DOC < prev   
Text File  |  1993-12-29  |  48KB  |  1,614 lines
  1.  
  2.  
  3.  
  4.  
  5.  
  6.   teste.doc - test.doc for lsodE
  7.  
  8.  
  9.   this file contains five demos and their outputs.  the file origin for
  10.   the first two files is /netlib/odepack/demo.Z at research.att.com and 
  11.   for imsl.for is page 647-648, math/library, imsl manual, version 1.1,
  12.   jan, 1989.  they are either trimmed or modified for lsode package.
  13.  
  14.  
  15.  
  16.  
  17.         contents
  18.         
  19.         1. demoe.c
  20.         2. demoe.for
  21.         3. output of demoe.c
  22.         4. output of demoe.for
  23.         5. ivpag.for
  24.         6. imsl.for
  25.         7. imsl.c
  26.         8. outputs of 5, 6, and 7
  27.         9. outputs from example.for and example.c
  28.            (example.c is in lsode.h)
  29.            (example.for is in lsode.c)
  30.  
  31.  
  32.  
  33.  
  34.  
  35. **************************************************************************
  36. **************************************************************************
  37. ********************************* demo.c *********************************
  38. **************************************************************************
  39. **************************************************************************
  40.  
  41.  
  42. /* -------------------------------------------------------------------- */
  43. /*  demoe.c was translated from the file "demo" that was downloaded     */
  44. /*  from netlib at research.att.com to test clsode.                     */
  45. /* -------------------------------------------------------------------- */
  46. /* demonstration program for the lsode package.                         */
  47. /*  this is the version of 13 august, 1981.                             */
  48. /*                                                                      */
  49. /*  this version is in double precision.                                */
  50. /*                                                                      */
  51. /*  the package is used to solve two simple problems,                   */
  52. /*  one with a full jacobian, the other with a banded jacobian,         */
  53. /*  with all 8 of the appropriate values of mf in each case.            */
  54. /*  if the errors are too large, or other difficulty occurs,            */
  55. /*  a warning message is printed.  all output is on unit lout = 6.      */
  56. /* -------------------------------------------------------------------- */
  57.  
  58. #include <stdlib.h>
  59. #include <math.h>
  60. #include "lsode.h"
  61.  
  62. #define   pd(a,b)  pd[nr*(b-1)+a]
  63.  
  64. void fcn1(int[], double, double[], double[]);
  65. void jac1(int[], double, double[], int, int, double[], int);
  66. void fcn2(int[], double, double[], double[]);
  67. void jac2(int[], double, double[], int, int, double[], int);
  68. double edit2(double[], double);
  69.  
  70.  
  71. int main()
  72. {
  73.   int    i, iopar, iopt, iout, istate, itask, itol, *iwork,
  74.          leniw, lenrw, liw, lrw, mband,  mf, miter, miter1, meth,
  75.          ml, mu, neq[2], nerr, nfe, nfea, nje, nout, nqu, nst;
  76.   double atole[2], dtout, er, erm, ero, hu, rtol[2], *rwork, t, tout,
  77.          tout1, *y;
  78.  
  79.  
  80.    tout1 = 1.39283880203e0;
  81.    dtout = 2.214773875e0;
  82.  
  83.    neq[1] = 2;
  84.    nerr = 0;
  85.    itol = 1;
  86.    rtol[1] = 0.0e0;
  87.    atole[1] = 1.0e-6;
  88.    lrw = 697;
  89.    liw = 45;
  90.    iwork = ivector(liw);
  91.    rwork = dvector(lrw);
  92.    y = dvector(neq[1]);
  93.    iopt = 0;
  94. /*                                                                      */
  95. /*  first problem ----------------------------------------------------- */
  96. /*                                                                      */
  97.    nout = 4;
  98.    printf("\f\n\n demonstration program for clsode package\n\n\n\n");
  99.    printf("\n problem 1..   van der pol oscillator..");
  100.    printf("\n  xdotdot - 3*(1 - x**2)*xdot + x = 0,");
  101.    printf("  x(0) = 2,  xdot(0) = 0");
  102.    printf("\n neq = %2d",neq[1]);
  103.    printf("\n itol = %3d   rtol = %10.1e  atol = %10.1e",
  104.              itol, rtol[1], atole[1]);
  105.  
  106.    for(meth=1;meth<=2;meth++) {
  107.        for(miter1=1;miter1<=4;miter1++) {
  108.             miter = miter1 - 1;
  109.             mf = 10*meth + miter;
  110.             printf("\n\n\n\n mf = %3d\n\n\n",mf);
  111.             printf("     t\t\t      x\t\t       xdot");
  112.             printf("          nq     h\n");
  113.             t = 0.0e0;
  114.             y[1] = 2.0e0;
  115.             y[2] = 0.0e0;
  116.             itask = 1;
  117.             istate = 1;
  118.             tout = tout1;
  119.             ero = 0.0e0;
  120.             for(iout=1;iout<=nout;iout++) {
  121.                  lsode(fcn1,neq,y,&t,tout,itol,rtol,atole,itask,&istate,
  122.                           iopt,rwork,lrw,iwork,liw,jac1,mf,stdout);
  123.                  hu = rwork[11];
  124.                  nqu = iwork[14];
  125.                  printf("\n %15.5e %16.5e %14.3e %5d %14.3e",
  126.                                  t, y[1], y[2], nqu, hu);
  127.                  if (istate < 0) break;
  128.                  iopar = iout - 2*(iout/2);
  129.                  if (iopar == 0) {
  130.                        er = fabs(y[1])/atole[1];
  131.                        ero = max(ero,er);
  132.                        if (er >= 1000.0e0) {
  133.                             printf("\n\n warning.. error ");
  134.                             printf("exceeds 1000 * tolerance\n\n");
  135.                             nerr++;
  136.                        }
  137.                  }
  138.                  tout += dtout;
  139.             }
  140.             if (istate < 0) nerr++;
  141.             nst = iwork[11];
  142.             nfe = iwork[12];
  143.             nje = iwork[13];
  144.             lenrw = iwork[17];
  145.             leniw = iwork[18];
  146.             nfea = nfe;
  147.             if (miter == 2) nfea = nfe - neq[1]*nje;
  148.             if (miter == 3) nfea = nfe - nje;
  149.             printf("\n\n final statistics for this run..");
  150.             printf("\n rwork size = %4d   iwork size = %4d",lenrw,leniw);
  151.             printf("\n number of steps = %5d",nst);
  152.             printf("\n number of f-s   = %5d",nfe);
  153.             printf("\n (excluding j-s) = %5d",nfea);
  154.             printf("\n number of j-s   = %5d",nje);
  155.             printf("\n error overrun =  %10.2e",ero);
  156.        }
  157.    }
  158. /*                                                                      */
  159. /*  second problem ---------------------------------------------------- */
  160. /*                                                                      */
  161.    neq[1] = 25;
  162.    y = (double *) realloc(y, (unsigned)(neq[1]+1)*sizeof(double));
  163.  
  164.    ml = 5;
  165.    mu = 0;
  166.    iwork[1] = ml;
  167.    iwork[2] = mu;
  168.    mband = ml + mu + 1;
  169.    nout = 5;
  170.    printf("\f\n\n demonstration program for clsode package\n\n\n\n");
  171.    printf("\n problem 2.. ydot = a * y , where,");
  172.    printf(" a is a banded lower triangular matrix");
  173.    printf("\n  derived from 2-d advection pde");
  174.    printf("\n neq = %3d   ml = %2d   mu = %2d",neq[1],ml,mu);
  175.    printf("\n itol = %3d   rtol = %10.1e   atol = %10.1e",
  176.                itol, rtol[1], atole[1]);
  177.  
  178.    for(meth=1;meth<=2;meth++) {
  179.        for(miter1=1;miter1<=6;miter1++) {
  180.             miter = miter1 - 1;
  181.             if (miter == 1 || miter == 2) goto skip_step;
  182.             mf = 10*meth + miter;
  183.             printf("\n\n\n\n mf = %3d\n\n\n",mf);
  184.             printf("     t\t\t   max.err          nq     h\n");
  185.             t = 0.0e0;
  186.             for(i=2;i<=neq[1];i++) y[i] = 0.0e0;
  187.             y[1] = 1.0e0;
  188.             itask = 1;
  189.             istate = 1;
  190.             tout = 0.01e0;
  191.             ero = 0.0e0;
  192.             for(iout=1;iout<=nout;iout++) {
  193.                  lsode(fcn2,neq,y,&t,tout,itol,rtol,atole,itask,&istate,
  194.                           iopt,rwork,lrw,iwork,liw,jac2,mf,stdout);
  195.                  erm = edit2(y,t);
  196.                  hu = rwork[11];
  197.                  nqu = iwork[14];
  198.                  printf("\n %15.5e %14.3e %5d %14.3e",t,erm,nqu,hu);
  199.                  if (istate < 0) break;
  200.                  er = erm/atole[1];
  201.                  ero = max(ero,er);
  202.                  if (er > 1000.0e0) {
  203.                        printf("\n\n warning.. error ");
  204.                        printf("exceeds 1000 * tolerance\n\n");
  205.                        nerr++;
  206.                  }
  207.                  tout *= 10.0e0;
  208.             }
  209.             if (istate < 0) nerr++;
  210.             nst = iwork[11];
  211.             nfe = iwork[12];
  212.             nje = iwork[13];
  213.             lenrw = iwork[17];
  214.             leniw = iwork[18];
  215.             nfea = nfe;
  216.             if (miter == 5) nfea = nfe - mband*nje;
  217.             if (miter == 3) nfea = nfe - nje;
  218.             printf("\n\n final statistics for this run..");
  219.             printf("\n rwork size = %4d   iwork size = %4d",lenrw,leniw);
  220.             printf("\n number of steps = %5d",nst);
  221.             printf("\n number of f-s   = %5d",nfe);
  222.             printf("\n (excluding j-s) = %5d",nfea);
  223.             printf("\n number of j-s   = %5d",nje);
  224.             printf("\n error overrun =  %10.2e",ero);
  225.   skip_step:  continue;
  226.       }
  227.    }
  228.    printf("\n\n\n\n number of errors encountered = %3d",nerr);
  229.    free(iwork);
  230.    free(rwork);
  231.    free(y);
  232.    return 1;
  233. }
  234.  
  235.  
  236. /*                                                                      */
  237. /*  derivative and jacobi matrix routine ------------------------------ */
  238. /*                                                                      */
  239.  
  240. void fcn1(int neq[], double t, double y[], double ydot[])
  241. {
  242.      ydot[1] = y[2];
  243.      ydot[2] = 3.0e0*(1.0e0 - y[1]*y[1])*y[2] - y[1];
  244.      return;
  245. }
  246.  
  247.  
  248. void jac1(int neq[], double t, double y[],
  249.                     int ml, int mu, double pd[], int nr)
  250. {
  251.       pd(1,1) = 0.0e0;
  252.       pd(1,2) = 1.0e0;
  253.       pd(2,1) = -6.0e0*y[1]*y[2] - 1.0e0;
  254.       pd(2,2) = 3.0e0*(1.0e0 - y[1]*y[1]);
  255.       return;
  256. }
  257.  
  258.  
  259. void fcn2(int neq[], double t, double y[], double ydot[])
  260. {
  261.       int i, j, k, ng= 5;
  262.       double alph1 = 1.0e0, alph2 = 1.0e0, d;
  263.  
  264.       for(j = 1;j<=ng;j++) {
  265.             for(i=1;i<=ng;i++) {
  266.                  k = i + (j - 1)*ng;
  267.                  d = -2.0e0*y[k];
  268.                  if (i != 1) d = d + y[k-1]*alph1;
  269.                  if (j != 1) d = d + y[k-ng]*alph2;
  270.                  ydot[k] = d;
  271.          }
  272.       }
  273.       return;
  274. }
  275.  
  276.  
  277. void jac2(int neq[], double t, double y[],
  278.                     int ml, int mu, double pd[], int nr)
  279. {
  280.       int j, mband, mu1, mu2, ng=5;
  281.       double alph1 = 1.0e0, alph2 = 1.0e0;
  282.  
  283.       mband = ml + mu + 1;
  284.       mu1 = mu + 1;
  285.       mu2 = mu + 2;
  286.       for(j=1;j<=neq[1];j++) {
  287.              pd(mu1,j) = -2.0e0;
  288.              pd(mu2,j) = alph1;
  289.              pd(mband,j) = alph2;
  290.       }
  291.       for(j=ng;j<=neq[1];j+=ng)  pd(mu2,j) = 0.0e0;
  292.       return;
  293. }
  294.  
  295.  
  296. double edit2 (double y[], double t)
  297. {
  298.     int i, j, k, ng=5;
  299.     double  alph1=1.0e0, alph2=1.0e0;
  300.     double  erm, a1, a2, er, ex, yt;
  301.  
  302.      erm = 0.0e0;
  303.      if (t == 0.0e0) return erm;
  304.      ex = 0.0e0;
  305.      if (t <= 30.0e0) ex = exp(-2.0e0*t);
  306.      a2 = 1.0e0;
  307.      for(j=1;j<=ng;j++) {
  308.          a1 = 1.0e0;
  309.          for(i=1;i<=ng;i++) {
  310.               k = i + (j - 1)*ng;
  311.               yt = pow(t,(i+j-2))*ex*a1*a2;
  312.               er = fabs(y[k]-yt);
  313.               erm = max(erm,er);
  314.               a1 *= alph1/(double)i;
  315.          }
  316.          a2 *= alph2/(double)j;
  317.     }
  318.     return erm;
  319. }
  320.  
  321. /* ---------------------------- end of demoe.c ------------------------ */
  322.  
  323.  
  324. **************************************************************************
  325. **************************************************************************
  326. ******************************* demoe.for ********************************
  327. **************************************************************************
  328. **************************************************************************
  329.  
  330. c-----------------------------------------------------------------------
  331. c demonstration program for the lsode package.
  332. c this is the version of 13 august, 1981.
  333. c
  334. c this version is in double precision.
  335. c
  336. c for computer systems requiring a program card, the following (with
  337. c the c in column 1 removed) may be used..
  338. c     program lsdem(lsout,tape6=lsout)
  339. c
  340. c the package is used to solve two simple problems,
  341. c one with a full jacobian, the other with a banded jacobian,
  342. c with all 8 of the appropriate values of mf in each case.
  343. c if the errors are too large, or other difficulty occurs,
  344. c a warning message is printed.  all output is on unit lout = 6.
  345. c-----------------------------------------------------------------------
  346.       external f1, jac1, f2, jac2
  347.       integer i, iopar, iopt, iout, istate, itask, itol, iwork,
  348.      1   leniw, lenrw, liw, lout, lrw, mband, meth, mf, miter, miter1,
  349.      2   ml, mu, neq, nerr, nfe, nfea, nje, nout, nqu, nst
  350.       double precision atol, dtout, er, erm, ero, hu, rtol, rwork, t,
  351.      1   tout, tout1, y
  352.       dimension y(25), rwork(697), iwork(45)
  353.       data lout/6/, tout1/1.39283880203d0/, dtout/2.214773875d0/
  354. c
  355.       nerr = 0
  356.       itol = 1
  357.       rtol = 0.0d0
  358.       atol = 1.0d-6
  359.       lrw = 697
  360.       liw = 45
  361.       iopt = 0
  362. c
  363. c first problem
  364. c
  365.       neq = 2
  366.       nout = 4
  367.       write (lout,110) neq,itol,rtol,atol
  368.  110  format(1h1/1x,40h demonstration program for lsode package////
  369.      1  1x,39h problem 1..   van der pol oscillator../
  370.      2  1x,38h  xdotdot - 3*(1 - x**2)*xdot + x = 0, ,
  371.      3  24h   x(0) = 2, xdot(0) = 0/
  372.      4  1x,6h neq =,i2/
  373.      5  1x,7h itol =,i3,9h   rtol =,e10.1,9h   atol =,e10.1//)
  374. c
  375.       do 195 meth = 1,2
  376.       do 190 miter1 = 1,4
  377.       miter = miter1 - 1
  378.       mf = 10*meth + miter
  379.       write (lout,120) mf
  380.  120  format(////1x,5h mf =,i3///
  381.      1  6x,1ht,15x,1hx,15x,4hxdot,7x,2hnq,6x,1hh/)
  382.       t = 0.0d0
  383.       y(1) = 2.0d0
  384.       y(2) = 0.0d0
  385.       itask = 1
  386.       istate = 1
  387.       tout = tout1
  388.       ero = 0.0d0
  389.       do 170 iout = 1,nout
  390.         call lsode(f1,neq,y,t,tout,itol,rtol,atol,itask,istate,
  391.      1     iopt,rwork,lrw,iwork,liw,jac1,mf)
  392.         hu = rwork(11)
  393.         nqu = iwork(14)
  394.         write (lout,140) t,y(1),y(2),nqu,hu
  395.  140    format(1x,e15.5,e16.5,e14.3,i5,e14.3)
  396.         if (istate .lt. 0) go to 175
  397.         iopar = iout - 2*(iout/2)
  398.         if (iopar .ne. 0) go to 170
  399.         er = dabs(y(1))/atol
  400.         ero = dmax1(ero,er)
  401.         if (er .lt. 1000.0d0) go to 170
  402.         write (lout,150)
  403.  150    format(//1x,41h warning.. error exceeds 1000 * tolerance//)
  404.         nerr = nerr + 1
  405.  170    tout = tout + dtout
  406.  175  continue
  407.       if (istate .lt. 0) nerr = nerr + 1
  408.       nst = iwork(11)
  409.       nfe = iwork(12)
  410.       nje = iwork(13)
  411.       lenrw = iwork(17)
  412.       leniw = iwork(18)
  413.       nfea = nfe
  414.       if (miter .eq. 2) nfea = nfe - neq*nje
  415.       if (miter .eq. 3) nfea = nfe - nje
  416.       write (lout,180) lenrw,leniw,nst,nfe,nfea,nje,ero
  417.  180  format(/1x,32h final statistics for this run../
  418.      1  1x,13h rwork size =,i4,15h   iwork size =,i4/
  419.      2  1x,18h number of steps =,i5/
  420.      3  1x,18h number of f-s   =,i5/
  421.      4  1x,18h (excluding j-s) =,i5/
  422.      5  1x,18h number of j-s   =,i5/
  423.      6  1x,16h error overrun =,e10.2)
  424.  190  continue
  425.  195  continue
  426. c
  427. c second problem
  428. c
  429.       neq = 25
  430.       ml = 5
  431.       mu = 0
  432.       iwork(1) = ml
  433.       iwork(2) = mu
  434.       mband = ml + mu + 1
  435.       nout = 5
  436.       write (lout,210) neq,ml,mu,itol,rtol,atol
  437.  210  format(1h1/1x,40h demonstration program for lsode package////
  438.      1  1x,33h problem 2.. ydot = a * y , where,
  439.      2  39h  a is a banded lower triangular matrix/
  440.      3  1x,32h  derived from 2-d advection pde/
  441.      4  1x,6h neq =,i3,7h   ml =,i2,7h   mu =,i2/
  442.      5  1x,7h itol =,i3,9h   rtol =,e10.1,9h   atol =,e10.1//)
  443.       do 295 meth = 1,2
  444.       do 290 miter1 = 1,6
  445.       miter = miter1 - 1
  446.       if (miter .eq. 1 .or. miter .eq. 2) go to 290
  447.       mf = 10*meth + miter
  448.       write (lout,220) mf
  449.  220  format(////1x,5h mf =,i3///
  450.      1  6x,1ht,13x,8hmax.err.,5x,2hnq,6x,1hh/)
  451.       t = 0.0d0
  452.       do 230 i = 2,neq
  453.  230    y(i) = 0.0d0
  454.       y(1) = 1.0d0
  455.       itask = 1
  456.       istate = 1
  457.       tout = 0.01d0
  458.       ero = 0.0d0
  459.       do 270 iout = 1,nout
  460.         call lsode(f2,neq,y,t,tout,itol,rtol,atol,itask,istate,
  461.      1     iopt,rwork,lrw,iwork,liw,jac2,mf)
  462.         call edit2(y,t,erm)
  463.         hu = rwork(11)
  464.         nqu = iwork(14)
  465.         write (lout,240) t,erm,nqu,hu
  466.  240    format(1x,e15.5,e14.3,i5,e14.3)
  467.         if (istate .lt. 0) go to 275
  468.         er = erm/atol
  469.         ero = dmax1(ero,er)
  470.         if (er .le. 1000.0d0) go to 270
  471.         write (lout,150)
  472.         nerr = nerr + 1
  473.  270    tout = tout*10.0d0
  474.  275  continue
  475.       if (istate .lt. 0) nerr = nerr + 1
  476.       nst = iwork(11)
  477.       nfe = iwork(12)
  478.       nje = iwork(13)
  479.       lenrw = iwork(17)
  480.       leniw = iwork(18)
  481.       nfea = nfe
  482.       if (miter .eq. 5) nfea = nfe - mband*nje
  483.       if (miter .eq. 3) nfea = nfe - nje
  484.       write (lout,180) lenrw,leniw,nst,nfe,nfea,nje,ero
  485.  290  continue
  486.  295  continue
  487.       write (lout,300) nerr
  488.  300  format(////1x,31h number of errors encountered =,i3)
  489.       stop
  490.       end
  491. c
  492.       subroutine f1 (neq, t, y, ydot)
  493.       integer neq
  494.       double precision t, y, ydot
  495.       dimension y(2), ydot(2)
  496.       ydot(1) = y(2)
  497.       ydot(2) = 3.0d0*(1.0d0 - y(1)*y(1))*y(2) - y(1)
  498.       return
  499.       end
  500.       subroutine jac1 (neq, t, y, ml, mu, pd, nrowpd)
  501.       integer neq, ml, mu, nrowpd
  502.       double precision t, y, pd
  503.       dimension y(2), pd(nrowpd,2)
  504.       pd(1,1) = 0.0d0
  505.       pd(1,2) = 1.0d0
  506.       pd(2,1) = -6.0d0*y(1)*y(2) - 1.0d0
  507.       pd(2,2) = 3.0d0*(1.0d0 - y(1)*y(1))
  508.       return
  509.       end
  510. c
  511.       subroutine f2 (neq, t, y, ydot)
  512.       integer neq, i, j, k, ng
  513.       double precision t, y, ydot, alph1, alph2, d
  514.       dimension y(1), ydot(1)
  515.       data alph1/1.0d0/, alph2/1.0d0/, ng/5/
  516.       do 10 j = 1,ng
  517.       do 10 i = 1,ng
  518.         k = i + (j - 1)*ng
  519.         d = -2.0d0*y(k)
  520.         if (i .ne. 1) d = d + y(k-1)*alph1
  521.         if (j .ne. 1) d = d + y(k-ng)*alph2
  522.  10     ydot(k) = d
  523.       return
  524.       end
  525. c      
  526.       subroutine jac2 (neq, t, y, ml, mu, pd, nrowpd)
  527.       integer neq, ml, mu, nrowpd, j, mband, mu1, mu2, ng
  528.       double precision t, y, pd, alph1, alph2
  529.       dimension y(1), pd(nrowpd,1)
  530.       data alph1/1.0d0/, alph2/1.0d0/, ng/5/
  531.       mband = ml + mu + 1
  532.       mu1 = mu + 1
  533.       mu2 = mu + 2
  534.       do 10 j = 1,neq
  535.         pd(mu1,j) = -2.0d0
  536.         pd(mu2,j) = alph1
  537.  10     pd(mband,j) = alph2
  538.       do 20 j = ng,neq,ng
  539.  20     pd(mu2,j) = 0.0d0
  540.       return
  541.       end
  542. c      
  543.       subroutine edit2 (y, t, erm)
  544.       integer i, j, k, ng
  545.       double precision y, t, erm, alph1, alph2, a1, a2, er, ex, yt
  546.       dimension y(25)
  547.       data alph1/1.0d0/, alph2/1.0d0/, ng/5/
  548.       erm = 0.0d0
  549.       if (t .eq. 0.0d0) return
  550.       ex = 0.0d0
  551.       if (t .le. 30.0d0) ex = dexp(-2.0d0*t)
  552.       a2 = 1.0d0
  553.       do 60 j = 1,ng
  554.         a1 = 1.0d0
  555.         do 50 i = 1,ng
  556.           k = i + (j - 1)*ng
  557.           yt = t**(i+j-2)*ex*a1*a2
  558.           er = dabs(y(k)-yt)
  559.           erm = dmax1(erm,er)
  560.           a1 = a1*alph1/dfloat(i)
  561.  50       continue
  562.         a2 = a2*alph2/dfloat(j)
  563.  60     continue
  564.       return
  565.       end
  566.       
  567. c ------------------------------ end of demo.for --------------------- 
  568.   
  569.  
  570. **************************************************************************
  571. **************************************************************************
  572. *************************** output of demoe.c ****************************
  573. **************************************************************************
  574. **************************************************************************
  575.  
  576.  demonstration program for clsode package
  577.  
  578.  
  579.  
  580.  
  581.  problem 1..   van der pol oscillator..
  582.   xdotdot - 3*(1 - x**2)*xdot + x = 0,  x(0) = 2,  xdot(0) = 0
  583.  neq =  2
  584.  itol =   1   rtol =    0.0e+00  atol =    1.0e-06
  585.  
  586.  
  587.  
  588.  mf =  10
  589.  
  590.  
  591.      t                x                xdot          nq     h
  592.  
  593.      1.39284e+00      1.68010e+00     -2.911e-01     3      1.229e-01
  594.      3.60761e+00     -7.79864e-05     -3.169e+00     5      2.171e-02
  595.      5.82239e+00     -1.68009e+00      2.911e-01     3      4.753e-02
  596.      8.03716e+00      1.16694e-04      3.169e+00     5      2.342e-02
  597.  
  598.  final statistics for this run..
  599.  rwork size =   52   iwork size =   20
  600.  number of steps =   297
  601.  number of f-s   =   352
  602.  (excluding j-s) =   352
  603.  number of j-s   =     0
  604.  error overrun =    1.17e+02
  605.  
  606.  
  607.  
  608.  mf =  11
  609.  
  610.  
  611.      t                x                xdot          nq     h
  612.  
  613.      1.39284e+00      1.68010e+00     -2.911e-01     5      1.205e-01
  614.      3.60761e+00     -1.77321e-05     -3.169e+00     5      1.872e-02
  615.      5.82239e+00     -1.68010e+00      2.911e-01     6      9.633e-02
  616.      8.03716e+00      2.58940e-05      3.169e+00     5      1.899e-02
  617.  
  618.  final statistics for this run..
  619.  rwork size =   58   iwork size =   22
  620.  number of steps =   203
  621.  number of f-s   =   281
  622.  (excluding j-s) =   281
  623.  number of j-s   =    29
  624.  error overrun =    2.59e+01
  625.  
  626.  
  627.  
  628.  mf =  12
  629.  
  630.  
  631.      t                x                xdot          nq     h
  632.  
  633.      1.39284e+00      1.68010e+00     -2.911e-01     5      1.205e-01
  634.      3.60761e+00     -1.77321e-05     -3.169e+00     5      1.872e-02
  635.      5.82239e+00     -1.68010e+00      2.911e-01     6      9.633e-02
  636.      8.03716e+00      2.58939e-05      3.169e+00     5      1.899e-02
  637.  
  638.  final statistics for this run..
  639.  rwork size =   58   iwork size =   22
  640.  number of steps =   203
  641.  number of f-s   =   339
  642.  (excluding j-s) =   281
  643.  number of j-s   =    29
  644.  error overrun =    2.59e+01
  645.  
  646.  
  647.  
  648.  mf =  13
  649.  
  650.  
  651.      t                x                xdot          nq     h
  652.  
  653.      1.39284e+00      1.68010e+00     -2.911e-01     5      7.386e-02
  654.      3.60761e+00      3.44007e-05     -3.169e+00     6      2.605e-02
  655.      5.82239e+00     -1.68011e+00      2.911e-01     4      1.328e-01
  656.      8.03716e+00     -5.91397e-05      3.169e+00     5      2.045e-02
  657.  
  658.  final statistics for this run..
  659.  rwork size =   56   iwork size =   20
  660.  number of steps =   198
  661.  number of f-s   =   315
  662.  (excluding j-s) =   289
  663.  number of j-s   =    26
  664.  error overrun =    5.91e+01
  665.  
  666.  
  667.  
  668.  mf =  20
  669.  
  670.      
  671.      t                x                xdot          nq     h
  672.  
  673.      1.39284e+00      1.68010e+00     -2.911e-01     5      5.490e-02
  674.      3.60761e+00     -5.65788e-05     -3.169e+00     5      1.429e-02
  675.      5.82239e+00     -1.68010e+00      2.911e-01     4      5.829e-02
  676.      8.03716e+00      1.03869e-04      3.169e+00     5      1.488e-02
  677.  
  678.  final statistics for this run..
  679.  rwork size =   38   iwork size =   20
  680.  number of steps =   289
  681.  number of f-s   =   321
  682.  (excluding j-s) =   321
  683.  number of j-s   =     0
  684.  error overrun =    1.04e+02
  685.  
  686.  
  687.  
  688.  mf =  21
  689.  
  690.  
  691.      t                x                xdot          nq     h
  692.  
  693.      1.39284e+00      1.68010e+00     -2.911e-01     5      6.762e-02
  694.      3.60761e+00     -4.89765e-05     -3.169e+00     5      1.413e-02
  695.      5.82239e+00     -1.68010e+00      2.911e-01     5      1.256e-01
  696.      8.03716e+00      9.68674e-05      3.169e+00     5      1.420e-02
  697.  
  698.  final statistics for this run..
  699.  rwork size =   44   iwork size =   22
  700.  number of steps =   262
  701.  number of f-s   =   345
  702.  (excluding j-s) =   345
  703.  number of j-s   =    30
  704.  error overrun =    9.69e+01
  705.  
  706.  
  707.  
  708.  mf =  22
  709.  
  710.  
  711.      t                x                xdot          nq     h
  712.  
  713.      1.39284e+00      1.68010e+00     -2.911e-01     5      6.762e-02
  714.      3.60761e+00     -4.89765e-05     -3.169e+00     5      1.413e-02
  715.      5.82239e+00     -1.68010e+00      2.911e-01     5      1.256e-01
  716.      8.03716e+00      9.68674e-05      3.169e+00     5      1.420e-02
  717.  
  718.  final statistics for this run..
  719.  rwork size =   44   iwork size =   22
  720.  number of steps =   262
  721.  number of f-s   =   405
  722.  (excluding j-s) =   345
  723.  number of j-s   =    30
  724.  error overrun =    9.69e+01
  725.  
  726.  
  727.  
  728.  mf =  23
  729.  
  730.  
  731.      t                x                xdot          nq     h
  732.  
  733.      1.39284e+00      1.68010e+00     -2.911e-01     5      7.089e-02
  734.      3.60761e+00     -4.67048e-05     -3.169e+00     5      1.390e-02
  735.      5.82239e+00     -1.68010e+00      2.911e-01     3      7.192e-02
  736.      8.03716e+00      5.47001e-05      3.169e+00     5      1.538e-02
  737.  
  738.  final statistics for this run..
  739.  rwork size =   42   iwork size =   20
  740.  number of steps =   271
  741.  number of f-s   =   414
  742.  (excluding j-s) =   383
  743.  number of j-s   =    31
  744.  error overrun =    5.47e+01
  745.  
  746.  
  747.  
  748.  demonstration program for clsode package
  749.  
  750.  
  751.  
  752.  
  753.  problem 2.. ydot = a * y , where, a is a banded lower triangular matrix
  754.   derived from 2-d advection pde
  755.  neq =  25   ml =  5   mu =  0
  756.  itol =   1   rtol =    0.0e+00   atol =    1.0e-06
  757.  
  758.  
  759.  
  760.  mf =  10
  761.  
  762.  
  763.      t             max.err          nq     h
  764.  
  765.      1.00000e-02      5.563e-07     2      7.660e-03
  766.      1.00000e-01      6.550e-06     3      2.494e-02
  767.      1.00000e+00      2.744e-06     4      5.198e-02
  768.      1.00000e+01      1.143e-06     3      1.169e-01
  769.      1.00000e+02      2.215e-06     2      2.615e-01
  770.  
  771.  final statistics for this run..
  772.  rwork size =  420   iwork size = 20
  773.  number of steps =   524
  774.  number of f-s   =   552
  775.  (excluding j-s) =   552
  776.  number of j-s   =     0
  777.  error overrun =    6.55e+00
  778.  
  779.  
  780.  
  781.  mf =  13
  782.  
  783.  
  784.      t             max.err          nq     h
  785.  
  786.      1.00000e-02      8.390e-07     2      9.493e-03
  787.      1.00000e-01      2.075e-06     3      2.497e-02
  788.      1.00000e+00      1.268e-04     3      1.677e-02
  789.      1.00000e+01      1.130e-05     3      3.853e-01
  790.      1.00000e+02      1.475e-06     2      1.486e+01
  791.  
  792.  final statistics for this run..
  793.  rwork size =  447   iwork size = 20
  794.  number of steps =   129
  795.  number of f-s   =   235
  796.  (excluding j-s) =   201
  797.  number of j-s   =    34
  798.  error overrun =    1.27e+02
  799.  
  800.  
  801.  
  802.  mf =  14
  803.  
  804.  
  805.      t             max.err          nq     h
  806.  
  807.      1.00000e-02      8.772e-07     2      9.648e-03
  808.      1.00000e-01      2.064e-06     3      2.497e-02
  809.      1.00000e+00      1.263e-06     5      9.348e-02
  810.      1.00000e+01      3.108e-07     6      4.423e-01
  811.      1.00000e+02      1.594e-08     2      2.909e+01
  812.  
  813.  final statistics for this run..
  814.  rwork size =  697   iwork size = 45
  815.  number of steps =    92
  816.  number of f-s   =   113
  817.  (excluding j-s) =   113
  818.  number of j-s   =    18
  819.  error overrun =    2.06e+00
  820.  
  821.  
  822.  
  823.  mf =  15
  824.  
  825.  
  826.      t             max.err          nq     h
  827.  
  828.      1.00000e-02      8.771e-07     2      9.648e-03
  829.      1.00000e-01      2.064e-06     3      2.497e-02
  830.      1.00000e+00      1.263e-06     5      9.348e-02
  831.      1.00000e+01      3.108e-07     6      4.423e-01
  832.      1.00000e+02      1.595e-08     2      2.909e+01
  833.  
  834.  final statistics for this run..
  835.  rwork size =  697   iwork size = 45
  836.  number of steps =    92
  837.  number of f-s   =   221
  838.  (excluding j-s) =   113
  839.  number of j-s   =    18
  840.  error overrun =    2.06e+00
  841.  
  842.  
  843.  
  844.  mf =  20
  845.  
  846.  
  847.      t             max.err          nq     h
  848.  
  849.      1.00000e-02      4.648e-07     2      4.827e-03
  850.      1.00000e-01      1.307e-06     3      1.477e-02
  851.      1.00000e+00      4.267e-06     5      6.347e-02
  852.      1.00000e+01      1.921e-06     4      3.514e-01
  853.      1.00000e+02      9.884e-08     1      4.526e-01
  854.  
  855.  final statistics for this run..
  856.  rwork size =  245   iwork size = 20
  857.  number of steps =   330
  858.  number of f-s   =   530
  859.  (excluding j-s) =   530
  860.  number of j-s   =     0
  861.  error overrun =    4.27e+00
  862.  
  863.  
  864.  
  865.  mf =  23
  866.  
  867.  
  868.      t             max.err          nq     h
  869.  
  870.      1.00000e-02      1.006e-06     2      5.983e-03
  871.      1.00000e-01      4.456e-07     3      1.459e-02
  872.      1.00000e+00      1.527e-06     5      7.378e-02
  873.      1.00000e+01      5.778e-07     4      3.235e-01
  874.      1.00000e+02      5.610e-12     1      4.747e+01
  875.  
  876.  final statistics for this run..
  877.  rwork size =  272   iwork size = 20
  878.  number of steps =   189
  879.  number of f-s   =   343
  880.  (excluding j-s) =   290
  881.  number of j-s   =    53
  882.  error overrun =    1.53e+00
  883.  
  884.  
  885.  
  886.  mf =  24
  887.  
  888.  
  889.      t             max.err          nq     h
  890.  
  891.      1.00000e-02      1.039e-06     2      6.081e-03
  892.      1.00000e-01      4.632e-07     3      1.459e-02
  893.      1.00000e+00      2.475e-06     5      6.658e-02
  894.      1.00000e+01      8.280e-07     5      3.912e-01
  895.      1.00000e+02      3.843e-10     1      1.084e+02
  896.  
  897.  final statistics for this run..
  898.  rwork size =  522   iwork size = 45
  899.  number of steps =   118
  900.  number of f-s   =   136
  901.  (excluding j-s) =   136
  902.  number of j-s   =    18
  903.  error overrun =    2.47e+00
  904.  
  905.  
  906.  
  907.  mf =  25
  908.  
  909.  
  910.      t             max.err          nq     h
  911.  
  912.      1.00000e-02      1.039e-06     2      6.081e-03
  913.      1.00000e-01      4.632e-07     3      1.459e-02
  914.      1.00000e+00      2.475e-06     5      6.658e-02
  915.      1.00000e+01      8.280e-07     5      3.912e-01
  916.      1.00000e+02      3.843e-10     1      1.084e+02
  917.  
  918.  final statistics for this run..
  919.  rwork size =  522   iwork size = 45
  920.  number of steps =   118
  921.  number of f-s   =   244
  922.  (excluding j-s) =   136
  923.  number of j-s   =    18
  924.  error overrun =    2.47e+00
  925.  
  926.  
  927.  
  928.  number of errors encountered =   0
  929.  
  930.  
  931. **************************************************************************
  932. **************************************************************************
  933. ************************** output of demoe.for ***************************
  934. **************************************************************************
  935. **************************************************************************
  936.  
  937.  
  938.  demonstration program for lsode package
  939.  
  940.  
  941.  
  942.  problem 1..   van der pol oscillator..
  943.   xdotdot - 3*(1 - x**2)*xdot + x = 0,   x(0) = 2, xdot(0) = 0
  944.  neq = 2
  945.  itol =  1   rtol =    .0E+00   atol =    .1E-05
  946.  
  947.  
  948.  
  949.  
  950.  
  951.  
  952.  mf = 10
  953.  
  954.  
  955.      t               x               xdot       nq      h
  956.  
  957.      .13928E+01      .16801E+01     -.291E+00    3      .123E+00
  958.      .36076E+01     -.77986E-04     -.317E+01    5      .217E-01
  959.      .58224E+01     -.16801E+01      .291E+00    3      .475E-01
  960.      .80372E+01      .11669E-03      .317E+01    5      .234E-01
  961.  
  962.  final statistics for this run..
  963.  rwork size =  52   iwork size =  20
  964.  number of steps =  297
  965.  number of f-s   =  352
  966.  (excluding j-s) =  352
  967.  number of j-s   =    0
  968.  error overrun =   .12E+03
  969.  
  970.  
  971.  
  972.  
  973.  mf = 11
  974.  
  975.  
  976.      t               x               xdot       nq      h
  977.  
  978.      .13928E+01      .16801E+01     -.291E+00    5      .121E+00
  979.      .36076E+01     -.17732E-04     -.317E+01    5      .187E-01
  980.      .58224E+01     -.16801E+01      .291E+00    6      .963E-01
  981.      .80372E+01      .25894E-04      .317E+01    5      .190E-01
  982.  
  983.  final statistics for this run..
  984.  rwork size =  58   iwork size =  22
  985.  number of steps =  203
  986.  number of f-s   =  281
  987.  (excluding j-s) =  281
  988.  number of j-s   =   29
  989.  error overrun =   .26E+02
  990.  
  991.  
  992.  
  993.  
  994.  mf = 12
  995.  
  996.  
  997.      t               x               xdot       nq      h
  998.  
  999.      .13928E+01      .16801E+01     -.291E+00    5      .121E+00
  1000.      .36076E+01     -.17732E-04     -.317E+01    5      .187E-01
  1001.      .58224E+01     -.16801E+01      .291E+00    6      .963E-01
  1002.      .80372E+01      .25894E-04      .317E+01    5      .190E-01
  1003.  
  1004.  final statistics for this run..
  1005.  rwork size =  58   iwork size =  22
  1006.  number of steps =  203
  1007.  number of f-s   =  339
  1008.  (excluding j-s) =  281
  1009.  number of j-s   =   29
  1010.  error overrun =   .26E+02
  1011.  
  1012.  
  1013.  
  1014.  
  1015.  mf = 13
  1016.  
  1017.  
  1018.      t               x               xdot       nq      h
  1019.  
  1020.      .13928E+01      .16801E+01     -.291E+00    5      .739E-01
  1021.      .36076E+01      .34401E-04     -.317E+01    6      .260E-01
  1022.      .58224E+01     -.16801E+01      .291E+00    4      .133E+00
  1023.      .80372E+01     -.59140E-04      .317E+01    5      .205E-01
  1024.  
  1025.  final statistics for this run..
  1026.  rwork size =  56   iwork size =  20
  1027.  number of steps =  198
  1028.  number of f-s   =  315
  1029.  (excluding j-s) =  289
  1030.  number of j-s   =   26
  1031.  error overrun =   .59E+02
  1032.  
  1033.  
  1034.  
  1035.  
  1036.  mf = 20
  1037.  
  1038.  
  1039.      t               x               xdot       nq      h
  1040.  
  1041.      .13928E+01      .16801E+01     -.291E+00    5      .549E-01
  1042.      .36076E+01     -.56579E-04     -.317E+01    5      .143E-01
  1043.      .58224E+01     -.16801E+01      .291E+00    4      .583E-01
  1044.      .80372E+01      .10387E-03      .317E+01    5      .149E-01
  1045.  
  1046.  final statistics for this run..
  1047.  rwork size =  38   iwork size =  20
  1048.  number of steps =  289
  1049.  number of f-s   =  321
  1050.  (excluding j-s) =  321
  1051.  number of j-s   =    0
  1052.  error overrun =   .10E+03
  1053.  
  1054.  
  1055.  
  1056.  
  1057.  mf = 21
  1058.  
  1059.  
  1060.      t               x               xdot       nq      h
  1061.  
  1062.      .13928E+01      .16801E+01     -.291E+00    5      .676E-01
  1063.      .36076E+01     -.48977E-04     -.317E+01    5      .141E-01
  1064.      .58224E+01     -.16801E+01      .291E+00    5      .126E+00
  1065.      .80372E+01      .96867E-04      .317E+01    5      .142E-01
  1066.  
  1067.  final statistics for this run..
  1068.  rwork size =  44   iwork size =  22
  1069.  number of steps =  262
  1070.  number of f-s   =  345
  1071.  (excluding j-s) =  345
  1072.  number of j-s   =   30
  1073.  error overrun =   .97E+02
  1074.  
  1075.  
  1076.  
  1077.  
  1078.  mf = 22
  1079.  
  1080.  
  1081.      t               x               xdot       nq      h
  1082.  
  1083.      .13928E+01      .16801E+01     -.291E+00    5      .676E-01
  1084.      .36076E+01     -.48977E-04     -.317E+01    5      .141E-01
  1085.      .58224E+01     -.16801E+01      .291E+00    5      .126E+00
  1086.      .80372E+01      .96867E-04      .317E+01    5      .142E-01
  1087.  
  1088.  final statistics for this run..
  1089.  rwork size =  44   iwork size =  22
  1090.  number of steps =  262
  1091.  number of f-s   =  405
  1092.  (excluding j-s) =  345
  1093.  number of j-s   =   30
  1094.  error overrun =   .97E+02
  1095.  
  1096.  
  1097.  
  1098.  
  1099.  mf = 23
  1100.  
  1101.  
  1102.      t               x               xdot       nq      h
  1103.  
  1104.      .13928E+01      .16801E+01     -.291E+00    5      .709E-01
  1105.      .36076E+01     -.46705E-04     -.317E+01    5      .139E-01
  1106.      .58224E+01     -.16801E+01      .291E+00    3      .719E-01
  1107.      .80372E+01      .54700E-04      .317E+01    5      .154E-01
  1108.  
  1109.  final statistics for this run..
  1110.  rwork size =  42   iwork size =  20
  1111.  number of steps =  271
  1112.  number of f-s   =  414
  1113.  (excluding j-s) =  383
  1114.  number of j-s   =   31
  1115.  error overrun =   .55E+02
  1116.  
  1117.  demonstration program for lsode package
  1118.  
  1119.  
  1120.  
  1121.  problem 2.. ydot = a * y , where  a is a banded lower triangular matrix
  1122.   derived from 2-d advection pde
  1123.  neq = 25   ml = 5   mu = 0
  1124.  itol =  1   rtol =    .0E+00   atol =    .1E-05
  1125.  
  1126.  
  1127.  
  1128.  
  1129.  
  1130.  mf = 10
  1131.  
  1132.  
  1133.      t             max.err.     nq      h
  1134.  
  1135.      .10000E-01      .556E-06    2      .766E-02
  1136.      .10000E+00      .655E-05    3      .249E-01
  1137.      .10000E+01      .274E-05    4      .520E-01
  1138.      .10000E+02      .114E-05    3      .117E+00
  1139.      .10000E+03      .221E-05    2      .262E+00
  1140.  
  1141.  final statistics for this run..
  1142.  rwork size = 420   iwork size =  20
  1143.  number of steps =  524
  1144.  number of f-s   =  552
  1145.  (excluding j-s) =  552
  1146.  number of j-s   =    0
  1147.  error overrun =   .65E+01
  1148.  
  1149.  
  1150.  
  1151.  
  1152.  mf = 13
  1153.  
  1154.  
  1155.      t             max.err.     nq      h
  1156.  
  1157.      .10000E-01      .839E-06    2      .949E-02
  1158.      .10000E+00      .208E-05    3      .250E-01
  1159.      .10000E+01      .127E-03    3      .168E-01
  1160.      .10000E+02      .113E-04    3      .385E+00
  1161.      .10000E+03      .147E-05    2      .149E+02
  1162.  
  1163.  final statistics for this run..
  1164.  rwork size = 447   iwork size =  20
  1165.  number of steps =  129
  1166.  number of f-s   =  235
  1167.  (excluding j-s) =  201
  1168.  number of j-s   =   34
  1169.  error overrun =   .13E+03
  1170.  
  1171.  
  1172.  
  1173.  
  1174.  mf = 14
  1175.  
  1176.  
  1177.      t             max.err.     nq      h
  1178.  
  1179.      .10000E-01      .877E-06    2      .965E-02
  1180.      .10000E+00      .206E-05    3      .250E-01
  1181.      .10000E+01      .126E-05    5      .935E-01
  1182.      .10000E+02      .311E-06    6      .442E+00
  1183.      .10000E+03      .159E-07    2      .291E+02
  1184.  
  1185.  
  1186.  final statistics for this run..
  1187.  rwork size = 697   iwork size =  45
  1188.  number of steps =   92
  1189.  number of f-s   =  113
  1190.  (excluding j-s) =  113
  1191.  number of j-s   =   18
  1192.  error overrun =   .21E+01
  1193.  
  1194.  
  1195.  
  1196.  
  1197.  mf = 15
  1198.  
  1199.  
  1200.      t             max.err.     nq      h
  1201.  
  1202.      .10000E-01      .877E-06    2      .965E-02
  1203.      .10000E+00      .206E-05    3      .250E-01
  1204.      .10000E+01      .126E-05    5      .935E-01
  1205.      .10000E+02      .311E-06    6      .442E+00
  1206.      .10000E+03      .160E-07    2      .291E+02
  1207.  
  1208.  final statistics for this run..
  1209.  rwork size = 697   iwork size =  45
  1210.  number of steps =   92
  1211.  number of f-s   =  221
  1212.  (excluding j-s) =  113
  1213.  number of j-s   =   18
  1214.  error overrun =   .21E+01
  1215.  
  1216.  
  1217.  
  1218.  
  1219.  mf = 20
  1220.  
  1221.  
  1222.      t             max.err.     nq      h
  1223.  
  1224.      .10000E-01      .465E-06    2      .483E-02
  1225.      .10000E+00      .131E-05    3      .148E-01
  1226.      .10000E+01      .427E-05    5      .635E-01
  1227.      .10000E+02      .192E-05    4      .351E+00
  1228.      .10000E+03      .988E-07    1      .453E+00
  1229.  
  1230.  final statistics for this run..
  1231.  rwork size = 245   iwork size =  20
  1232.  number of steps =  330
  1233.  number of f-s   =  530
  1234.  (excluding j-s) =  530
  1235.  number of j-s   =    0
  1236.  error overrun =   .43E+01
  1237.  
  1238.  
  1239.  
  1240.  
  1241.  mf = 23
  1242.  
  1243.  
  1244.      t             max.err.     nq      h
  1245.  
  1246.      .10000E-01      .101E-05    2      .598E-02
  1247.      .10000E+00      .446E-06    3      .146E-01
  1248.      .10000E+01      .153E-05    5      .738E-01
  1249.      .10000E+02      .578E-06    4      .324E+00
  1250.      .10000E+03      .561E-11    1      .475E+02
  1251.  
  1252.  final statistics for this run..
  1253.  rwork size = 272   iwork size =  20
  1254.  number of steps =  189
  1255.  number of f-s   =  343
  1256.  (excluding j-s) =  290
  1257.  number of j-s   =   53
  1258.  error overrun =   .15E+01
  1259.  
  1260.  
  1261.  
  1262.  
  1263.  mf = 24
  1264.  
  1265.  
  1266.      t             max.err.     nq      h
  1267.  
  1268.      .10000E-01      .104E-05    2      .608E-02
  1269.      .10000E+00      .463E-06    3      .146E-01
  1270.      .10000E+01      .247E-05    5      .666E-01
  1271.      .10000E+02      .828E-06    5      .391E+00
  1272.      .10000E+03      .384E-09    1      .108E+03
  1273.  
  1274.  final statistics for this run..
  1275.  rwork size = 522   iwork size =  45
  1276.  number of steps =  118
  1277.  number of f-s   =  136
  1278.  (excluding j-s) =  136
  1279.  number of j-s   =   18
  1280.  error overrun =   .25E+01
  1281.  
  1282.  
  1283.  
  1284.  
  1285.  mf = 25
  1286.  
  1287.  
  1288.      t             max.err.     nq      h
  1289.  
  1290.      .10000E-01      .104E-05    2      .608E-02
  1291.      .10000E+00      .463E-06    3      .146E-01
  1292.      .10000E+01      .247E-05    5      .666E-01
  1293.      .10000E+02      .828E-06    5      .391E+00
  1294.      .10000E+03      .384E-09    1      .108E+03
  1295.  
  1296.  final statistics for this run..
  1297.  rwork size = 522   iwork size =  45
  1298.  number of steps =  118
  1299.  number of f-s   =  244
  1300.  (excluding j-s) =  136
  1301.  number of j-s   =   18
  1302.  error overrun =   .25E+01
  1303.  
  1304.  
  1305.  
  1306.  
  1307.  number of errors encountered =  0
  1308. Stop - Program terminated.
  1309.   
  1310.  
  1311. **************************************************************************
  1312. **************************************************************************
  1313. ******************************* ivpag.for ********************************
  1314. **************************************************************************
  1315. **************************************************************************
  1316.  
  1317.  
  1318. C  IVPAG.FOR
  1319. C  FILE ORIGIN:
  1320. C    IMSL USER'S MANUAL, IMSL MATH/LIBRARY, VERSION 1.1, JAN, 1989
  1321. C    PAGE 647 - 648
  1322. C
  1323.       INTEGER   NEQ, NPARAM
  1324.       PARAMETER (NEQ=3, NPARAM=50)
  1325. C
  1326.       INTEGER   IDO, IEND, IMETH, INORM, NOUT
  1327.       REAL      A(1,1), FCN, FCNJ, HINIT, PARAM(NPARAM), TOL, X,
  1328.      &          XEND, Y(NEQ)
  1329.       EXTERNAL  FCN, FCNJ, IVPAG, SSET, UMACH
  1330. C                                 Initialize
  1331.       HINIT = 1.0E-03
  1332.       INORM = 2
  1333.       IMETH = 1
  1334.       CALL SSET (NPARAM, 0.0, PARAM, 1)
  1335.       PARAM(1)  = HINIT
  1336.       PARAM(10) = INORM 
  1337.       PARAM(12) = IMETH
  1338. C
  1339.       IDO  = 1
  1340.       X    = 0.0
  1341.       Y(1) = 0.0
  1342.       Y(2) = 1.0
  1343.       Y(3) = 1.0
  1344.       TOL  = 1.0E-6
  1345. C                                 Write title
  1346.       CALL UMACH (2, NOUT)
  1347.       WRITE (NOUT,99998)
  1348. C                                 Integrate ODE
  1349.       DO 10 IEND=1, 10
  1350.            XEND = IEND
  1351.            CALL IVPAG (IDO, NEQ, FCN, FCNJ, A, X, XEND, TOL, PARAM, Y)
  1352.            WRITE (NOUT,99999) X, Y
  1353.    10 CONTINUE
  1354. C                                 Finish up
  1355.       IDO = 3
  1356.       CALL IVPAG (IDO, NEQ, FCN, FCNJ, A, X, XEND, TOL, PARAM, Y)
  1357. C
  1358. 99998 FORMAT (11X, 'X', 14X, 'Y(1)', 11X, 'Y(2)', 11X, 'Y(3)')
  1359. 99999 FORMAT (4F15.5)
  1360.       END
  1361. C
  1362.       SUBROUTINE FCN (NEQ, X, Y, YPRIME)
  1363.       INTEGER    NEQ
  1364.       DIMENSION  Y(NEQ), YPRIME(NEQ)
  1365. C
  1366.       YPRIME(1) = Y(2)*Y(3)
  1367.       YPRIME(2) = -Y(1)*Y(3)
  1368.       YPRIME(3) = -0.51D0*Y(1)*Y(2)
  1369.       RETURN
  1370.       END
  1371. C
  1372.       SUBROUTINE FCNJ(NEQ, X, Y, DYPDY)
  1373.       INTEGER    NEQ
  1374.       DIMENSION  X, Y(NEQ), DYPDY(*)
  1375. C                                 This subroutine is never called
  1376.       RETURN
  1377.       END
  1378.  
  1379.  
  1380. **************************************************************************
  1381. **************************************************************************
  1382. ********************************* imsl.for *******************************
  1383. **************************************************************************
  1384. **************************************************************************
  1385.  
  1386. c
  1387. c  imsl.for
  1388. c  file origin:
  1389. c    imsl user's manual, imsl math/library, version 1.1, jan, 1989
  1390. c    page 647 - 648
  1391. c
  1392.       double precision atol,rtol,rwork,t,tout,y
  1393.       integer i,iopt,istate,itask,itol,iwork,liw,lrw,mf,neq
  1394.       dimension iwork(23),rwork(68),y(3)
  1395.       external fcn,jac
  1396. c
  1397. c connect unit 6 for output to the screen.
  1398. c
  1399.       open (6,file='con')
  1400. c
  1401. c enable display of error messages on unit 6.
  1402. c
  1403.       call xsetf(1)
  1404.       call xsetun(6)
  1405. c
  1406. c initialization of lsode call list parameters.
  1407. c
  1408.       neq = 3
  1409.       y(1) = 0.0d0
  1410.       y(2) = 1.0d0
  1411.       y(3) = 1.0d0
  1412. c
  1413.       t = 0.d0
  1414.       tout = 0.0d0
  1415. c
  1416.       itol = 1
  1417.       rtol = 1.d-4
  1418.       atol = 1.d-06
  1419. c
  1420.       itask = 1
  1421.       istate = 1
  1422.       iopt = 1
  1423.       lrw = 68
  1424.       liw = 23
  1425. c
  1426.       mf = 10
  1427.       write(6,24)
  1428. c
  1429. c solution and output of results.
  1430. c
  1431.          do 10, i = 0,10
  1432.             call lsode(fcn,neq,y,t,tout,itol,rtol,atol,itask,istate,
  1433.      *                  iopt,rwork,lrw,iwork,liw,fcnj,mf)
  1434.             write(6,20) t,y
  1435.             if (istate .lt. 0) then
  1436.                 write (6,60) istate
  1437.                 write (6,30) iwork(11)
  1438.                 write (6,40) iwork(12)
  1439.                 stop
  1440.             endif
  1441.             tout = tout + 1.0d0
  1442.    10    continue
  1443. c      
  1444.       write (6,30) iwork(11)
  1445.       write (6,40) iwork(12)
  1446. c
  1447.    20 format (4f15.5)
  1448.    24 format (/,6x,'time',14x,'y(1)',16x,'y(2)',16x,'y(3)'/)
  1449.    25 format (3x,f9.5,6x,3(f12.7,6x))
  1450.    30 format (/,1x,'number of steps = ',i4)
  1451.    40 format (1x,'number of fcn evaluations =',i4)
  1452.    60 format (/,1x,'error halt:  istate = ',i2)
  1453. c
  1454.       end
  1455. c
  1456. c
  1457. c  fcn defines the system of odes.
  1458. c
  1459.       subroutine fcn(neq, x, y, yprime)
  1460. c
  1461.       double precision x, y, yprime
  1462.       integer neq
  1463.       dimension y(neq), yprime(neq)
  1464. c
  1465.       yprime(1) = y(2)*y(3)
  1466.       yprime(2) = -y(1)*y(3)
  1467.       yprime(3) = -0.51d0*y(1)*y(2)
  1468.       return
  1469.       end
  1470. c
  1471. c
  1472. c  jac computes the jacobian matrix.
  1473. c  this subroutine is never called.
  1474. c
  1475.       subroutine fcnj(neq,x,y,ml,mu,dypdy,nr)
  1476. c
  1477.       double precision x, y, dypdy
  1478.       integer ml,mu,neq,nr
  1479.       dimension y(1), dypdy(nr,1)
  1480. c
  1481.       return
  1482.       end
  1483.  
  1484.  
  1485.       
  1486. **************************************************************************
  1487. **************************************************************************
  1488. ********************************** imsl.c ********************************
  1489. **************************************************************************
  1490. **************************************************************************
  1491.  
  1492.  
  1493. /*  imsl.c                                                             */
  1494. /*  file origin:                                                       */
  1495. /*    imsl user's manual, imsl math/library, version 1.1, jan, 1989    */
  1496. /*    page 647 - 648                                                   */
  1497.  
  1498.  
  1499. #include <math.h>
  1500. #include <stdlib.h>
  1501. #include "lsode.h"
  1502.  
  1503.  
  1504. void fcn(int[], double, double[], double[]);
  1505.  
  1506. int main()
  1507. {
  1508.       double atole[1+1], rtol[1+1], *rwork,t,tout,*y;
  1509.       int i,iopt,istate,itask,itol,*iwork,liw,lrw,mf,neq[1+1];
  1510.  
  1511.       neq[1] = 3;
  1512.       y = dvector(neq[1]);
  1513.       y[1] = 0.0e0;
  1514.       y[2] = 1.0e0;
  1515.       y[3] = 1.0e0;
  1516.  
  1517.       t = 0.0e0;
  1518.       tout = 0.0e0;
  1519.  
  1520.       itol = 1;
  1521.       rtol[1] = 1.e-04;
  1522.       atole[1] = 1.e-06;
  1523.  
  1524.       itask = 1;
  1525.       istate = 1;
  1526.       iopt = 1;
  1527.       lrw = 68;
  1528.       liw = 23;
  1529.       rwork = dvector(lrw);
  1530.       iwork = ivector(liw);
  1531.  
  1532.       mf = 10;
  1533.  
  1534. /*                                                                     */
  1535. /*  integration step - call lsode. ----------------------------------- */
  1536. /*                                                                     */
  1537.         for(i=0; i<=10; i++) {
  1538.              lsode(fcn,neq,y,&t,tout,itol,rtol,atole,itask,&istate,
  1539.                          iopt,rwork,lrw,iwork,liw,NULL,mf,stdout);
  1540.              if (istate < 0) {
  1541.                    printf("\n error halt:  istate = %3d",istate);
  1542.                    break;
  1543.              }
  1544.              printf("\n%15.5f %15.5f %15.5f %15.5f",t,y[1],y[2],y[3]);
  1545.              tout += 1.0e0;
  1546.       }
  1547.  
  1548.       printf("\n\n  number of steps = %d",iwork[11]);
  1549.       printf("\n  number of fcn evaluations = %d",iwork[12]);
  1550.  
  1551.       free(iwork);
  1552.       free(rwork);
  1553.       free(y);
  1554.       return 0;
  1555. }
  1556.  
  1557.  
  1558. void  fcn(int neq[], double t, double y[], double yprime[])
  1559. {
  1560.      yprime[1] = y[2]*y[3];
  1561.      yprime[2] = -y[1]*y[3];
  1562.      yprime[3] = -0.51e0*y[1]*y[2];
  1563.      return;
  1564. }
  1565.  
  1566. **************************************************************************
  1567. **************************************************************************
  1568. **************** outputs of ivpag.for, imsl.for, and imsl.c **************
  1569. **************************************************************************
  1570. **************************************************************************
  1571.  
  1572.   OUTPUT OF IVPAG.FOR  (SINGLE PRECISION RESULT)
  1573.                   
  1574.            X                Y(1)             Y(2)             Y(3)
  1575.         1.00000          0.80220          0.59705          0.81963
  1576.         2.00000          0.99537         -0.09615          0.70336
  1577.         3.00000          0.64141         -0.76720          0.88892
  1578.         4.00000         -0.26961         -0.96296          0.98129
  1579.         5.00000         -0.91172         -0.41079          0.75899
  1580.         6.00000         -0.95750          0.28841          0.72967
  1581.         7.00000         -0.42877          0.90342          0.95197
  1582.         8.00000          0.51093          0.85963          0.93106
  1583.         9.00000          0.97567          0.21925          0.71730
  1584.        10.00000          0.87790         -0.47885          0.77906
  1585.        
  1586.   output of imsl.for (double precision result) n-s = 82   n-f = 89 
  1587.      
  1588.           time            y(1)              y(2)              y(3)
  1589.        1.00000          .80220            .59710            .81966   
  1590.        2.00000          .99528           -.09618            .70341   
  1591.        3.00000          .64130           -.76714            .88895   
  1592.        4.00000         -.26932           -.96319            .98143   
  1593.        5.00000         -.91187           -.41084            .75898   
  1594.        6.00000         -.95751            .28848            .72971   
  1595.        7.00000         -.42891            .90333            .95195   
  1596.        8.00000          .51058            .86030            .93137   
  1597.        9.00000          .97596            .21953            .71729   
  1598.       10.00000          .87823           -.47878            .77905   
  1599.   
  1600.   output of imsl.c  (double precision result)  n-s = 82   n-f = 89
  1601.         
  1602.            time             y[1]             y[2]             y[3]
  1603.         1.00000          0.80220          0.59710          0.81966
  1604.         2.00000          0.99528         -0.09618          0.70341
  1605.         3.00000          0.64130         -0.76714          0.88895
  1606.         4.00000         -0.26932         -0.96319          0.98143
  1607.         5.00000         -0.91187         -0.41084          0.75898
  1608.         6.00000         -0.95751          0.28848          0.72971
  1609.         7.00000         -0.42891          0.90333          0.95195
  1610.         8.00000          0.51058          0.86030          0.93137
  1611.         9.00000          0.97596          0.21953          0.71729
  1612.        10.00000          0.87823         -0.47878          0.77905
  1613.   
  1614.