home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1994 #1 / monster.zip / monster / PROG_C / C_LSODE.ZIP / EXAMPLE.C < prev    next >
Text File  |  1994-01-03  |  5KB  |  131 lines
  1.  
  2.  
  3. /*  example.c                                                            */
  4. /* --------------------------------------------------------------------- */
  5. /*  example problem.                                                     */
  6. /*     - it is the same program that f-lsode has in the statements       */
  7. /*                                                                       */
  8. /*  the following is a simple example problem, with the coding           */
  9. /*  needed for its solution by lsode.  the problem is from chemical      */
  10. /*  kinetics, and consists of the following three rate equations..       */
  11. /*      dy1/dt = -.04*y1 + 1.e4*y2*y3                                    */
  12. /*      dy2/dt = .04*y1 - 1.e4*y2*y3 - 3.e7*y2**2                        */
  13. /*      dy3/dt = 3.e7*y2**2                                              */
  14. /*  on the interval from t = 0.0 to t = 4.e10, with initial conditions   */
  15. /*  y1 = 1.0, y2 = y3 = 0.  the problem is stiff.                        */
  16. /*                                                                       */
  17. /*  the following coding solves this problem with lsode, using mf = 21   */
  18. /*  and printing results at t = .4, 4., ..., 4.e10.  it uses             */
  19. /*  itol = 2 and atol much smaller for y2 than y1 or y3 because          */
  20. /*  y2 has much smaller values.                                          */
  21. /*  at the end of the run, statistical quantities of interest are        */
  22. /*  printed (see optional outputs in the full description below).        */
  23. /* --------------------------------------------------------------------- */
  24.  
  25. #include <stdlib.h>
  26. #include "lsode.h"
  27.  
  28.  
  29.  
  30. void fcn(int[], double, double[], double[]);
  31. void jac(int[], double, double[], int, int, double[], int);
  32.  
  33.  
  34. int main()
  35. {
  36.     double atole[3+1],rtol[2],*rwork,t,tout,*y ;
  37.     int itol,itask,istate,iopt,lrw,liw,mf,*iwork,neq[2];
  38.     int i;
  39.  
  40.       itol = 2;
  41.       lrw = 58;
  42.       liw = 23;
  43.       neq[1] = 3;
  44.  
  45.       iwork = ivector(liw);
  46.       rwork = dvector(lrw);
  47.       y = dvector(neq[1]);
  48.  
  49.       y[1] = 1.0e0;
  50.       y[2] = 0.0e0;
  51.       y[3] = 0.0e0;
  52.       t = 0.0e0;
  53.       tout = 0.40e0;
  54.  
  55.       rtol[1] = 1.e-4;
  56.       atole[1] = 1.e-6;
  57.       atole[2] = 1.e-10;
  58.       atole[3] = 1.e-6;
  59.  
  60.       itask = 1;
  61.       istate = 1;
  62.       iopt = 0;
  63.       mf = 21;
  64.  
  65.       for(i=1;i<=12;i++) {
  66.             lsode(fcn,neq,y,&t,tout,itol,rtol,atole,itask,&istate,
  67.                       iopt,rwork,lrw,iwork,liw,jac,mf,stdout);
  68.             printf("\n at t =%11.4e  y =%14.6e  %14.6e  %14.6e",
  69.                        t, y[1], y[2], y[3]);
  70.             if (istate < 0) {
  71.                  printf("\n\n error halt.. istate =%3d",istate);
  72.                  break;
  73.             }
  74.             tout *= 10.0e0;
  75.       }
  76.  
  77.       printf("\n\n no. steps = %4d  no. f-s = %4d  no. j-s =%4d\n\n",
  78.                        iwork[11], iwork[12], iwork[13]);
  79.       free(iwork);
  80.       free(rwork);
  81.       free(y);
  82.       return 1;
  83. }
  84.  
  85.  
  86. void fcn(int neq[], double t, double y[], double ydot[])
  87. {
  88.       ydot[1] = -.04e0*y[1] + 1.e4*y[2]*y[3];
  89.       ydot[3] = 3.e7*y[2]*y[2];
  90.       ydot[2] = -ydot[1] - ydot[3];
  91.       return;
  92. }
  93.  
  94.  
  95.  
  96.  
  97. #define   pd(a,b)  pd[nr*(b-1)+a]
  98. void jac(int neq[],double t,double y[],int ml,int mu,double pd[],int nr)
  99. {
  100.       pd(1,1) = -0.04e0;                      /*  pd[1][1]  */
  101.       pd(1,2) = 1.e4*y[3];                    /*  pd[1][2]  */
  102.       pd(1,3) = 1.e4*y[2];                    /*  pd[1][3]  */
  103.       pd(2,1) = 0.04e0;                       /*  pd[2][1]  */
  104.       pd(2,3) = -pd(1,3);                     /*  pd[2][3]  */
  105.       pd(3,2) = 6.e7*y[2];                    /*  pd[3][2]  */
  106.       pd(2,2) = -pd(1,2) - pd(3,2);           /*  pd[2][2]  */
  107.       return;
  108. }
  109. /* --------------------------------------------------------------------- */
  110. /*  the output(*) of this program is as follows..                        */
  111. /*                                                                       */
  112. /*  at t = 4.0000e-01  y =  9.851726e-01   3.386406e-05   1.479357e-02   */
  113. /*  at t = 4.0000e+00  y =  9.055142e-01   2.240418e-05   9.446344e-02   */
  114. /*  at t = 4.0000e+01  y =  7.158050e-01   9.184616e-06   2.841858e-01   */
  115. /*  at t = 4.0000e+02  y =  4.504846e-01   3.222434e-06   5.495122e-01   */
  116. /*  at t = 4.0000e+03  y =  1.831701e-01   8.940379e-07   8.168290e-01   */
  117. /*  at t = 4.0000e+04  y =  3.897016e-02   1.621193e-07   9.610297e-01   */
  118. /*  at t = 4.0000e+05  y =  4.935213e-03   1.983756e-08   9.950648e-01   */
  119. /*  at t = 4.0000e+06  y =  5.159269e-04   2.064759e-09   9.994841e-01   */
  120. /*  at t = 4.0000e+07  y =  5.306413e-05   2.122677e-10   9.999469e-01   */
  121. /*  at t = 4.0000e+08  y =  5.494530e-06   2.197824e-11   9.999945e-01   */
  122. /*  at t = 4.0000e+09  y =  5.129458e-07   2.051784e-12   9.999995e-01   */
  123. /*  at t = 4.0000e+10  y = -7.170560e-08  -2.868224e-13   1.000000e+00   */
  124. /*                                                                       */
  125. /*  no. steps =  330  no. f-s =  405  no. j-s =  69                      */
  126. /*                                                                       */
  127. /*  *. the output is the same as that of f-lsode (see teste.doc).        */
  128. /*     however, it may vary slightly from compiler to compiler.          */
  129. /* --------------------------------------------------------------------- */
  130.  
  131.