home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / c_lsode.zip / IMSL.C < prev    next >
Text File  |  1993-09-15  |  2KB  |  77 lines
  1.  
  2. /*---------------------------------------------------------------------*/
  3. /*  imsl.c                                                             */
  4. /*  file origin:                                                       */
  5. /*    imsl user's manual, imsl math/library, version 1.1, jan, 1989    */
  6. /*    page 647 - 648                                                   */
  7. /*---------------------------------------------------------------------*/
  8.  
  9.  
  10. #include <math.h>
  11. #include <stdlib.h>
  12. #include "lsode.h"
  13.  
  14.  
  15. void fcn(int[], double, double[], double[]);
  16.  
  17. int main()
  18. {
  19.       double atole[1+1], rtol[1+1], *rwork,t,tout,*y;
  20.       int i,iopt,istate,itask,itol,*iwork,liw,lrw,mf,neq[1+1];
  21.  
  22.       neq[1] = 3;
  23.       y = dvector(neq[1]);
  24.       y[1] = 0.0e0;
  25.       y[2] = 1.0e0;
  26.       y[3] = 1.0e0;
  27.  
  28.       t = 0.0e0;
  29.       tout = 0.0e0;
  30.  
  31.       itol = 1;
  32.       rtol[1] = 1.e-04;
  33.       atole[1] = 1.e-06;
  34.  
  35.       itask = 1;
  36.       istate = 1;
  37.       iopt = 1;
  38.       lrw = 68;
  39.       liw = 23;
  40.       rwork = dvector(lrw);
  41.       iwork = ivector(liw);
  42.  
  43.       mf = 10;
  44.  
  45. /*  integration step - call lsode. ----------------------------------- */
  46.  
  47.         for(i=0; i<=10; i++) {
  48.              lsode(fcn,neq,y,&t,tout,itol,rtol,atole,itask,&istate,
  49.                          iopt,rwork,lrw,iwork,liw,NULL,mf,stdout);
  50.              if (istate < 0) {
  51.                    printf("\n error halt:  istate = %3d",istate);
  52.                    break;
  53.              }
  54.              printf("\n%15.5f  %15.5f  %15.5f  %15.5f",t,y[1],y[2],y[3]);
  55.              tout += 1.0e0;
  56.       }
  57.  
  58.       printf("\n\n  number of steps = %d",iwork[11]);
  59.       printf("\n  number of fcn evaluations = %d",iwork[12]);
  60.  
  61.       free(iwork);
  62.       free(rwork);
  63.       free(y);
  64.       return 0;
  65. }
  66.  
  67.  
  68.  
  69. void  fcn(int neq[], double t, double y[], double yprime[])
  70. {
  71.      yprime[1] = y[2]*y[3];
  72.      yprime[2] = -y[1]*y[3];
  73.      yprime[3] = -0.51e0*y[1]*y[2];
  74.      return;
  75. }
  76.  
  77.