home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / c_lsode.zip / DEMOE.C next >
Text File  |  1994-01-03  |  13KB  |  307 lines
  1.  
  2. /*  demo.c                                                               */
  3. /*  file origin:                                                         */
  4. /*    /netlib/odepack/demo.Z                                             */
  5. /*    node address: research.att.com                                     */
  6.  
  7. /* --------------------------------------------------------------------- */
  8. /*  see test.doc for both fortran and c sources and outputs.             */
  9. /* --------------------------------------------------------------------- */
  10. /*                                                                       */
  11. /* demonstration program for the lsode package.                          */
  12. /*  this is the version of 13 august, 1981.                              */
  13. /*                                                                       */
  14. /*  this version is in double precision.                                 */
  15. /*                                                                       */
  16. /*  the package is used to solve two simple problems,                    */
  17. /*  one with a full jacobian, the other with a banded jacobian,          */
  18. /*  with all 8 of the appropriate values of mf in each case.             */
  19. /*  if the errors are too large, or other difficulty occurs,             */
  20. /*  a warning message is printed.                                        */
  21. /*                                                                       */
  22. /*  this program has been tested with following compilers.               */
  23. /*    bc++ 3.1 for win/dos                                               */
  24. /*    djgpp gcc 1.10                                                     */
  25. /*    gcc/2 2.3.3 for os/2                                               */
  26. /*    ibm c set++ of os/2                                                */
  27. /*    msc/c++ 7.0 for win/dos                                            */
  28. /*    power c for dos                                                    */
  29. /*    c89 for dec ultrix                                                 */
  30. /*                                                                       */
  31. /*  about the result....                                                 */
  32. /*  + bc++ and ibm c++ generate exactly the same results with demo.c,    */
  33. /*    however they are slightly different with that of ms fortran 5.1.   */
  34. /*    ms fortran 5.1 produced two different results when with and/or     */ 
  35. /*    without time optimizations. you have to turn the optimization      */
  36. /*    off to obtain the same result you will have with c compilers.      */
  37. /*  + although you could have the same results from f77 and c89, they    */
  38. /*    are not identical to one from any dos or os/2 compiler.            */
  39. /*                                                                       */
  40. /*  recommendation for compiling....                                     */
  41. /*    ms fortran 5.1  - turn time optimization off.                      */
  42. /*    borland c++ 3.1 - no optimization for common subexpressions        */
  43. /*                                                                       */
  44. /* --------------------------------------------------------------------- */
  45.  
  46.  
  47. #include <stdlib.h>       
  48. #include <math.h>
  49. #include "lsode.h"
  50.  
  51.  
  52. void fcn1(int [], double, double [], double []);
  53. void jac1(int [], double, double [], int, int, double [], int);
  54. void fcn2(int [], double, double [], double []);
  55. void jac2(int [], double, double [], int, int, double [], int);
  56. double edit2(double [], double);
  57.  
  58. int main()
  59. {
  60.   int    i, iopar, iopt, iout, istate, itask, itol, *iwork,
  61.          leniw, lenrw, liw, lrw, mband,  mf, miter, miter1, meth,
  62.          ml, mu, neq[2], nerr, nfe, nfea, nje, nout, nqu, nst;
  63.   double atole[2], dtout, er, erm, ero, hu, rtol[2], *rwork, t, tout,
  64.          tout1, *y;
  65.  
  66.      tout1 = 1.39283880203e0;
  67.      dtout = 2.214773875e0;
  68.  
  69.      neq[1] = 2;
  70.      nerr = 0;
  71.      itol = 1;
  72.      rtol[1] = 0.0e0;
  73.      atole[1] = 1.0e-6;
  74.      lrw = 697;
  75.      liw = 45;
  76.      iwork = ivector(liw);
  77.      rwork = dvector(lrw);
  78.      y = dvector(neq[1]);
  79.      iopt = 0;
  80. /*                                                                       */
  81. /*  first problem ------------------------------------------------------ */
  82. /*                                                                       */
  83.      nout = 4;
  84.      printf("\f\n\n demonstration program for clsode package\n\n\n\n");
  85.      printf("\n problem 1..   van der pol oscillator..");
  86.      printf("\n  xdotdot - 3*(1 - x**2)*xdot + x = 0,");
  87.      printf("  x(0) = 2,  xdot(0) = 0");
  88.      printf("\n neq = %2d",neq[1]);
  89.      printf("\n itol = %3d   rtol = %10.1e  atol = %10.1e",
  90.                   itol, rtol[1], atole[1]);
  91.  
  92.      for(meth=1;meth<=2;meth++) {
  93.          for(miter1=1;miter1<=4;miter1++) {
  94.              miter = miter1 - 1;
  95.              mf = 10*meth + miter;
  96.              printf("\n\n\n\n mf = %3d\n\n\n",mf);
  97.              printf("     t\t\t      x\t\t       xdot");
  98.              printf("          nq     h\n");
  99.              t = 0.0e0;
  100.              y[1] = 2.0e0;
  101.              y[2] = 0.0e0;
  102.              itask = 1;
  103.              istate = 1;
  104.              tout = tout1;
  105.              ero = 0.0e0;
  106.              for(iout=1;iout<=nout;iout++) {
  107.                   lsode(fcn1,neq,y,&t,tout,itol,rtol,atole,itask,&istate,
  108.                           iopt,rwork,lrw,iwork,liw,jac1,mf,stdout);
  109.                   hu = rwork[11];
  110.                   nqu = iwork[14];
  111.                   printf("\n %15.5e %16.5e %14.3e %5d %14.3e",
  112.                                  t, y[1], y[2], nqu, hu);
  113.                   if (istate < 0) break;
  114.                   iopar = iout - 2*(iout/2);
  115.                   if (iopar == 0) {
  116.                         er = fabs(y[1])/atole[1];
  117.                         ero = max(ero,er);
  118.                         if (er >= 1000.0e0) {
  119.                              printf("\n\n\n warning.. error ");
  120.                              printf("exceeds 1000 * tolerance\n\n");
  121.                              nerr = nerr + 1;
  122.                         }
  123.                   }
  124.                   tout += dtout;
  125.              }
  126.              if (istate < 0) nerr++;
  127.              nst = iwork[11];
  128.              nfe = iwork[12];
  129.              nje = iwork[13];
  130.              lenrw = iwork[17];
  131.              leniw = iwork[18];
  132.              nfea = nfe;
  133.              if (miter == 2) nfea = nfe - neq[1]*nje;
  134.              if (miter == 3) nfea = nfe - nje;
  135.              printf("\n\n final statistics for this run..");
  136.              printf("\n rwork size = %4d   iwork size = %4d",lenrw,leniw);
  137.              printf("\n number of steps = %5d",nst);
  138.              printf("\n number of f-s   = %5d",nfe);
  139.              printf("\n (excluding j-s) = %5d",nfea);
  140.              printf("\n number of j-s   = %5d",nje);
  141.              printf("\n error overrun =  %10.2e",ero);
  142.          }
  143.      }
  144. /*                                                                       */
  145. /*  second problem ----------------------------------------------------- */
  146. /*                                                                       */
  147.      neq[1] = 25;
  148.      y = (double *) realloc(y, (unsigned)(neq[1]+1)*sizeof(double));
  149.      ml = 5;
  150.      mu = 0;
  151.      iwork[1] = ml;
  152.      iwork[2] = mu;
  153.      mband = ml + mu + 1;
  154.      nout = 5;
  155.      printf("\f\n\n demonstration program for clsode package\n\n\n");
  156.      printf("\n\n problem 2.. ydot = a * y , where,");
  157.      printf(" a is a banded lower triangular matrix");
  158.      printf("\n  derived from 2-d advection pde");
  159.      printf("\n neq = %3d   ml = %2d   mu = %2d",neq[1],ml,mu);
  160.      printf("\n itol = %3d   rtol = %10.1e   atol = %10.1e",
  161.                     itol, rtol[1], atole[1]);
  162.  
  163.     for(meth=1;meth<=2;meth++) {
  164.         for(miter1=1;miter1<=6;miter1++) {
  165.              miter = miter1 - 1;
  166.              if (miter != 1 && miter != 2) {
  167.                   mf = 10*meth + miter;
  168.                   printf("\n\n\n\n mf = %3d\n\n\n     ",mf);
  169.                   printf("t\t\t   max.err          nq     h\n");
  170.                   t = 0.0e0;
  171.                   for(i=2;i<=neq[1];i++) y[i] = 0.0e0;
  172.                   y[1] = 1.0e0;
  173.                   itask = 1;
  174.                   istate = 1;
  175.                   tout = 0.01e0;
  176.                   ero = 0.0e0;
  177.                   for(iout=1;iout<=nout;iout++) {
  178.                         lsode(fcn2,neq,y,&t,tout,itol,rtol,atole,
  179.                                 itask,&istate,iopt,rwork,lrw,
  180.                                  iwork,liw,jac2,mf,stdout);
  181.                         erm = edit2(y,t);
  182.                         hu = rwork[11];
  183.                         nqu = iwork[14];
  184.                         printf("\n %15.5e %14.3e %5d %14.3e",t,erm,nqu,hu);
  185.                         if (istate < 0) break;
  186.                         er = erm/atole[1];
  187.                         ero = max(ero,er);
  188.                         if (er > 1000.0e0) {
  189.                              printf("\n\n\n warning.. error ");
  190.                              printf("exceeds 1000 * tolerance\n\n");
  191.                              nerr++;
  192.                         }
  193.                         tout *= 10.0e0;
  194.                   }
  195.                   if (istate < 0) nerr++;
  196.                   nst = iwork[11];
  197.                   nfe = iwork[12];
  198.                   nje = iwork[13];
  199.                   lenrw = iwork[17];
  200.                   leniw = iwork[18];
  201.                   nfea = nfe;
  202.                   if (miter == 5) nfea = nfe - mband*nje;
  203.                   if (miter == 3) nfea = nfe - nje;
  204.                   printf("\n\n final statistics for this run..");
  205.                   printf("\n rwork size = %4d   iwork size = %d",
  206.                              lenrw,leniw);
  207.                   printf("\n number of steps = %5d",nst);
  208.                   printf("\n number of f-s   = %5d",nfe);
  209.                   printf("\n (excluding j-s) = %5d",nfea);
  210.                   printf("\n number of j-s   = %5d",nje);
  211.                   printf("\n error overrun =  %10.2e",ero);
  212.              }
  213.         }
  214.     }
  215.     printf("\n\n\n number of errors encountered = %3d",nerr);
  216.     free(iwork);
  217.     free(rwork);
  218.     free(y);
  219.     return 1;
  220. }
  221.  
  222.  
  223. /*                                                                       */
  224. /*  derivative and jacobi matrix routine ------------------------------- */
  225. /*                                                                       */
  226.  
  227. void fcn1(int neq[], double t, double y[], double ydot[])
  228. {
  229.       ydot[1] = y[2];
  230.       ydot[2] = 3.0e0*(1.0e0 - y[1]*y[1])*y[2] - y[1];
  231.       return;
  232. }
  233.  
  234. #define   pd(a,b)  pd[nr*(b-1)+a]
  235. void jac1(int neq[], double t, double y[],
  236.               int ml, int mu, double pd[], int nr)
  237. {
  238.       pd(1,1) = 0.0e0;
  239.       pd(1,2) = 1.0e0;
  240.       pd(2,1) = -6.0e0*y[1]*y[2] - 1.0e0;
  241.       pd(2,2) = 3.0e0*(1.0e0 - y[1]*y[1]);
  242.       return;
  243. }
  244.  
  245.  
  246. void fcn2(int neq[], double t, double y[], double ydot[])
  247. {
  248.     int i, j, k, ng= 5;
  249.     double alph1 = 1.0e0, alph2 = 1.0e0, d;
  250.  
  251.       for(j = 1;j<=ng;j++) {
  252.             for(i=1;i<=ng;i++) {
  253.                  k = i + (j - 1)*ng;
  254.                  d = -2.0e0*y[k];
  255.                  if (i != 1) d += y[k-1]*alph1;
  256.                  if (j != 1) d += y[k-ng]*alph2;
  257.                  ydot[k] = d;
  258.             }
  259.       }
  260.       return;
  261. }
  262.  
  263.  
  264. void jac2(int neq[], double t, double y[],
  265.                     int ml, int mu, double pd[], int nr)
  266. {
  267.     int j, mband, mu1, mu2, ng=5;
  268.     double alph1 = 1.0e0, alph2 = 1.0e0;
  269.  
  270.       mband = ml + mu + 1;
  271.       mu1 = mu + 1;
  272.       mu2 = mu + 2;
  273.       for(j=1;j<=neq[1];j++) {
  274.             pd(mu1,j) = -2.0e0;
  275.             pd(mu2,j) = alph1;
  276.             pd(mband,j) = alph2;
  277.       }
  278.       for(j=ng;j<=neq[1];j+=ng) pd(mu2,j) = 0.0e0;
  279.       return;
  280. }
  281.  
  282.  
  283. double edit2 (double y[], double t)
  284. {
  285.     int  i, j, k, ng=5;
  286.     double  alph1=1.0e0, alph2=1.0e0;
  287.     double  a1, a2, er, ex, yt, erm;
  288.  
  289.       erm = 0.0e0;
  290.       if (t == 0.0e0) return erm;
  291.       ex = 0.0e0;
  292.       if (t <= 30.0e0) ex = exp(-2.0e0*t);
  293.       a2 = 1.0e0;
  294.       for(j=1;j<=ng;j++) {
  295.            a1 = 1.0e0;
  296.            for(i=1;i<=ng;i++) {
  297.                  k = i + (j - 1)*ng;
  298.                  yt = pow(t,(i+j-2))*ex*a1*a2;
  299.                  er = fabs(y[k]-yt);
  300.                  erm = max(erm,er);
  301.                  a1 *= alph1/(double)i;
  302.            }
  303.            a2 *= alph2/(double)j;
  304.       }
  305.       return erm;
  306. }
  307.