home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 7 / 07.iso / c / c129 / 1.ddi / INTEGRAT.C < prev    next >
Encoding:
C/C++ Source or Header  |  1988-01-28  |  4.7 KB  |  217 lines
  1. # include <math.h>
  2.  
  3.  
  4. float IFunction(float y)
  5. {
  6.    float _fret;
  7.  
  8.    _fret = exp(y);
  9.    return(_fret);
  10. }
  11.  
  12.  
  13. float YPrime(float x,float y)
  14. {
  15.    float _fret;
  16.  
  17.    _fret = exp(2*x)+y;
  18.    return(_fret);
  19. }
  20.  
  21.  
  22. void IntegratFunction(int numseg,float xl,float xh,float *ans)
  23. {
  24.    int i;
  25.    int k;
  26.    int segcntr;
  27.    float strtpnt;
  28.    float endpnt;
  29.    float area1;
  30.    float area2;
  31.    float area;
  32.    char even;
  33.    float segwidth;
  34.    segwidth = (xh - xl) / numseg;
  35.    endpnt = xh;
  36.    area1 = 0.0;
  37.    area2 = 0.0;
  38.    area = 0.0;
  39.    if ( fmod(numseg,2) != 0 ) {
  40.       area1 = 3.0 / 8.0 * segwidth * (IFunction(endpnt - 3.0 * segwidth) +
  41.               3.0 * IFunction(endpnt - 2.0 * segwidth) +
  42.               3.0 * IFunction(endpnt - segwidth) + IFunction(endpnt));
  43.       endpnt = endpnt - 3.0 * segwidth;
  44.    }
  45.    else {
  46.       area1 = 0.0;
  47.    }
  48.    if ( numseg != 3 ) {
  49.       strtpnt = xl;
  50.       while ( (strtpnt < endpnt - segwidth) ) {
  51.          area2 = area2 + 1.0 / 3.0 * (segwidth * (IFunction(strtpnt) + 4.0 *
  52.             IFunction(strtpnt + segwidth) + IFunction(strtpnt + 2.0 * segwidth)));
  53.          strtpnt = strtpnt + 2.0 * segwidth;
  54.       }
  55.    }
  56.    area = area1 + area2;
  57.    (*ans) = area;
  58. }
  59.  
  60. void IntegratVector(float y[],float samper,int xl,int xh,float *ans)
  61. {
  62.    int numseg;
  63.    int i;
  64.    int strtpnt;
  65.    int endpnt;
  66.    int segcntr;
  67.    float area1;
  68.    float area2;
  69.    float area;
  70.    char even;
  71.    float segwidth;
  72.  
  73.    numseg = xh - xl;
  74.    segwidth = samper;
  75.    endpnt = xh;
  76.    area1 = 0.0;
  77.    area2 = 0.0;
  78.    area = 0.0;
  79.    if ( fmod(numseg,2) != 0 ) {
  80.       area1 = 3.0 / 8.0 * segwidth * (y[endpnt - 3] + 3.0 *
  81.               y[endpnt - 2] + 3.0 * y[endpnt - 1] + y[endpnt]);
  82.       endpnt = endpnt - 3;
  83.    }
  84.    else {
  85.       area1 = 0.0;
  86.    }
  87.    if ( numseg != 3 ) {
  88.       strtpnt = xl;
  89.       do {
  90.          area2 = area2 + 1.0 / 3.0 * segwidth * (y[strtpnt] + 4.0
  91.                  * y[strtpnt + 1] + y[strtpnt + 2]);
  92.          strtpnt = strtpnt + 2;
  93.       } while ( ! (strtpnt == (endpnt)) );
  94.    }
  95.    area = area1 + area2;
  96.    (*ans) = area;
  97. }
  98.  
  99. void RungeKutta(float XKnown,float YKnown,float x,float et,float h,
  100.                 float HMin,float *YEst,char *status)
  101. {
  102.    # define NearZero 1.0e-20
  103.    float xc2;
  104.    float xc3;
  105.    float xc4;
  106.    float xc5;
  107.    float xc6;
  108.    float dc21;
  109.    float dc31;
  110.    float dc32;
  111.    float dc41;
  112.    float dc42;
  113.    float dc43;
  114.    float dc51;
  115.    float dc52;
  116.    float dc53;
  117.    float dc54;
  118.    float dc61;
  119.    float dc62;
  120.    float dc63;
  121.    float dc64;
  122.    float dc65;
  123.    float yc1;
  124.    float yc3;
  125.    float yc4;
  126.    float yc5;
  127.    float yc6;
  128.    float ec1;
  129.    float ec3;
  130.    float ec4;
  131.    float ec5;
  132.    float ec6;
  133.    float d1;
  134.    float d2;
  135.    float d3;
  136.    float d4;
  137.    float d5;
  138.    float d6;
  139.    float error;
  140.    float XEst;
  141.    float EtDiv32;
  142.    int st;
  143.  
  144.    xc2 = 0.25;
  145.    xc3 = 3.0 / 8.0;
  146.    xc4 = 12.0 / 13.0;
  147.    xc5 = 1.0;
  148.    xc6 = 0.5;
  149.    dc21 = 0.25;
  150.    dc31 = 3.0 / 32.0;
  151.    dc32 = 9.0 / 32.0;
  152.    dc41 = 1932.0 / 2197.0;
  153.    dc42 = 7200.0 / 2197.0;
  154.    dc43 = 7296.0 / 2197.0;
  155.    dc51 = 439.0 / 216.0;
  156.    dc52 = 8.0;
  157.    dc53 = 3680.0 / 513.0;
  158.    dc54 = 845.0 / 4104.0;
  159.    dc61 = 8.0 / 27.0;
  160.    dc62 = 2.0;
  161.    dc63 = 3544.0 / 2565.0;
  162.    dc64 = 1859.0 / 4104.0;
  163.    dc65 = 11.0 / 40.0;
  164.    yc1 = 25.0 / 216.0;
  165.    yc3 = 1408.0 / 2565.0;
  166.    yc4 = 2197.0 / 4104.0;
  167.    yc5 = 0.20;
  168.    ec1 = 1.0 / 360.0;
  169.    ec3 = 128.0 / 4275.0;
  170.    ec4 = 2197.0 / 75240.0;
  171.    ec5 = 0.02;
  172.    ec6 = 2.0 / 55.0;
  173.    EtDiv32 = et / 32.0;
  174.    XEst = XKnown;
  175.    (*YEst) = YKnown;
  176.    st = 0;
  177.    while ( st == 0 ) {
  178.       d1 = YPrime(XEst,(*YEst));
  179.       d2 = YPrime(XEst + h * xc2,(*YEst) + h * dc21 * d1);
  180.       d3 = YPrime(XEst + h * xc3,(*YEst) + h * (dc31 * d1 + dc32 * d2));
  181.       d4 = YPrime(XEst + h * xc4,(*YEst) + h * (dc41 * d1 - dc42 * d2 + dc43 * d3));
  182.       d5 = YPrime(XEst + h * xc5,(*YEst) + h * (dc51 * d1 - dc52 * d2 + dc53 * d3 - dc54 * d4));
  183.       d6 = YPrime(XEst + h * xc6,(*YEst) + h * (-dc61 * d1 + dc62 * d2 - dc63 * d3 + dc64 * d4 - dc65 * d5));
  184.       error = h * (ec1 * d1 - ec3 * d3 - ec4 * d4 + ec5 * d5 + ec6 * d6);
  185.  
  186.       if ( error < (h * et) ) {
  187.  
  188.          XEst = XEst + h;
  189.          (*YEst) = (*YEst) + h * (yc1 * d1 + yc3 * d3 + yc4 * d4 - yc5 * d5);
  190.          if ( fabs(x - XEst) > NearZero ) {
  191.             if ( error <= h * EtDiv32 ) {
  192.                h = h * 2.0;
  193.             }
  194.             if ( h > (x - XEst) ) {
  195.                h = x - XEst;
  196.             }
  197.          }
  198.          else {
  199.             st = 1;
  200.          }
  201.       }
  202.       else {
  203.          h = h / 2;
  204.          if ( h < HMin ) {
  205.             st = 2;
  206.          }
  207.       }
  208.       if ( st == 2 ) {
  209.          (*status) = 0;
  210.       }
  211.       else {
  212.          (*status) = 1;
  213.       }
  214.    }
  215. }
  216.  
  217.