program

This is a discussion on program within the Windows Programming forums, part of the Platform Specific Boards category; How do I write a program to solve the integral (e^x)*(x^4)/(e^x-1) from a to b where a is the lower ...

  1. #1
    Registered User
    Join Date
    Aug 2014
    Posts
    4

    program

    How do I write a program to solve the integral (e^x)*(x^4)/(e^x-1) from a to b where a is the lower limit and b is the upper limit input by the user. I am having serious problems writing the code and I need the program urgently.

  2. #2
    and the hat of wrongness Salem's Avatar
    Join Date
    Aug 2001
    Location
    The edge of the known universe
    Posts
    32,553
    Well the first thing you should do is post your current attempt.
    If you dance barefoot on the broken glass of undefined behaviour, you've got to expect the occasional cut.
    If at first you don't succeed, try writing your phone number on the exam paper.
    I support http://www.ukip.org/ as the first necessary step to a free Europe.

  3. #3
    Registered User
    Join Date
    Jun 2005
    Posts
    6,307
    What do you mean by "solve the integral"? As you've worded it, your question has no meaning.

    If you mean "integrate your expression (the integrand) with respect to x, over the interval [a,b]" then, personally, I'd find a way to simply the integrand symbolically, for example using Maple. Depending on the integrand (I haven't examined yours closely) you may be able to compute the integral analytically - and you're better doing that if at all possible. Failing that, simplifying the integrand may give you an expression in terms of elementary or special functions (in the sense of Abramowitz and Stegun's Handbook of Mathematical Functions), since there are good quality libraries (some even available free) to compute many of those functions.

    If it's not possible analytically or with the help of special functions (or if the integral is expressed in terms of a special function that is not in your chosen library) you can fall back on an algorithm for numerical integration. The word you want to google for is "quadrature", which literally means "a method or algorithm for numerically approximating a definite integral". There are plenty of types of quadrature rules, including Newton-Cotes formulae, the trapezoidal rule (which is simple and often effective, so a lot of people use it, even when it doesn't work well), Simpson's rule, and others.

    It's a bit hard (at least on a casual look) to say what algorithm would work effectively in your case - your integrand has a division by zero when x is zero, which will tend to upset quadrature techniques if a and b are of different signs (one positive, one negative).

    If you want more specific help, show what you've tried.
    Right 98% of the time, and don't care about the other 3%.

  4. #4
    Registered User
    Join Date
    Aug 2014
    Posts
    4
    well here is what I tried. I first integrated this through a math software and now I have this : 4x^3 Li2(e^x) - 12x^2 Li3(e^x) + 24x Li4(e^x) - 24 Li5(e^x) + x^4 log(1-e^x) + constant.
    here is then what I tried :
    Code:
    /*                            polylog.c
     *
     *    Polylogarithms
     *
     *
     *
     * SYNOPSIS:
     */
      double x, y, polylog();
      int n;
     
      y = polylog( n, x );
     
     
     /* The polylogarithm of order n is defined by the series
     
     *
     *              inf   k
     *               -   x
     *  Li (x)  =    >   ---  .
     *    n          -     n
     *              k=1   k
     *
     *
     *  For x = 1,
     *
     *               inf
     *                -    1
     *   Li (1)  =    >   ---   =  Riemann zeta function (n)  .
     *     n          -     n
     *               k=1   k
     *
     *
     *  When n = 2, the function is the dilogarithm, related to Spence's integral:
     *
     *                 x                      1-x
     *                 -                        -
     *                | |  -ln(1-t)            | |  ln t
     *   Li (x)  =    |    -------- dt    =    |    ------ dt    =   spence(1-x) .
     *     2        | |       t              | |    1 - t
     *               -                        -
     *                0                        1
     *
     *
     *  See also the program cpolylog.c for the complex polylogarithm,
     *  whose definition is extended to x > 1.
     *
     *  
     *
     * ACCURACY:
     *
     *                      Relative error:
     * arithmetic   domain   n   # trials      peak         rms
     *    IEEE      0, 1     2     50000      6.2e-16     8.0e-17
     *    IEEE      0, 1     3    100000      2.5e-16     6.6e-17
     *    IEEE      0, 1     4     30000      1.7e-16     4.9e-17
     *    IEEE      0, 1     5     30000      5.1e-16     7.8e-17
     *
     */
    
    
    
    
    
    #include "cmath"
    extern double PI;
    
    
    /* polylog(4, 1-x) = zeta(4) - x zeta(3) + x^2 A4(x)/B4(x)
       0 <= x <= 0.125
       Theoretical peak absolute error 4.5e-18  */
    #if UNK
    static double A4[13] = {
     3.056144922089490701751E-2,
     3.243086484162581557457E-1,
     2.877847281461875922565E-1,
     7.091267785886180663385E-2,
     6.466460072456621248630E-3,
     2.450233019296542883275E-4,
     4.031655364627704957049E-6,
     2.884169163909467997099E-8,
     8.680067002466594858347E-11,
     1.025983405866370985438E-13,
     4.233468313538272640380E-17,
     4.959422035066206902317E-21,
     1.059365867585275714599E-25,
    };
    static double B4[12] = {
      /* 1.000000000000000000000E0, */
     2.821262403600310974875E0,
     1.780221124881327022033E0,
     3.778888211867875721773E-1,
     3.193887040074337940323E-2,
     1.161252418498096498304E-3,
     1.867362374829870620091E-5,
     1.319022779715294371091E-7,
     3.942755256555603046095E-10,
     4.644326968986396928092E-13,
     1.913336021014307074861E-16,
     2.240041814626069927477E-20,
     4.784036597230791011855E-25,
    };
    #endif
    #if DEC
    static short A4[52] = {
    0036772,0056001,0016601,0164507,
    0037646,0005710,0076603,0176456,
    0037623,0054205,0013532,0026476,
    0037221,0035252,0101064,0065407,
    0036323,0162231,0042033,0107244,
    0035200,0073170,0106141,0136543,
    0033607,0043647,0163672,0055340,
    0031767,0137614,0173376,0072313,
    0027676,0160156,0161276,0034203,
    0025347,0003752,0123106,0064266,
    0022503,0035770,0160173,0177501,
    0017273,0056226,0033704,0132530,
    0013403,0022244,0175205,0052161,
    };
    static short B4[48] = {
      /*0040200,0000000,0000000,0000000, */
    0040464,0107620,0027471,0071672,
    0040343,0157111,0025601,0137255,
    0037701,0075244,0140412,0160220,
    0037002,0151125,0036572,0057163,
    0035630,0032452,0050727,0161653,
    0034234,0122515,0034323,0172615,
    0032415,0120405,0123660,0003160,
    0030330,0140530,0161045,0150177,
    0026002,0134747,0014542,0002510,
    0023134,0113666,0035730,0035732,
    0017723,0110343,0041217,0007764,
    0014024,0007412,0175575,0160230,
    };
    #endif
    #if IBMPC
    static short A4[52] = {
    0x3d29,0x23b0,0x4b80,0x3f9f,
    0x7fa6,0x0fb0,0xc179,0x3fd4,
    0x45a8,0xa2eb,0x6b10,0x3fd2,
    0x8d61,0x5046,0x2755,0x3fb2,
    0x71d4,0x2883,0x7c93,0x3f7a,
    0x37ac,0x118c,0x0ecf,0x3f30,
    0x4b5c,0xfcf7,0xe8f4,0x3ed0,
    0xce99,0x9edf,0xf7f1,0x3e5e,
    0xc710,0xdc57,0xdc0d,0x3dd7,
    0xcd17,0x54c8,0xe0fd,0x3d3c,
    0x7fe8,0x1c0f,0x677f,0x3c88,
    0x96ab,0xc6f8,0x6b92,0x3bb7,
    0xaa8e,0x9f50,0x6494,0x3ac0,
    };
    static short B4[48] = {
      /*0x0000,0x0000,0x0000,0x3ff0,*/
    0x2e77,0x05e7,0x91f2,0x4006,
    0x37d6,0x2570,0x7bc9,0x3ffc,
    0x5c12,0x9821,0x2f54,0x3fd8,
    0x4bce,0xa7af,0x5a4a,0x3fa0,
    0xfc75,0x4a3a,0x06a5,0x3f53,
    0x7eb2,0xa71a,0x94a9,0x3ef3,
    0x00ce,0xb4f6,0xb420,0x3e81,
    0xba10,0x1c44,0x182b,0x3dfb,
    0x40a9,0xe32c,0x573c,0x3d60,
    0x077b,0xc77b,0x92f6,0x3cab,
    0xe1fe,0x6851,0x721c,0x3bda,
    0xbc13,0x5f6f,0x81e1,0x3ae2,
    };
    #endif
    #if MIEEE
    static short A4[52] = {
    0x3f9f,0x4b80,0x23b0,0x3d29,
    0x3fd4,0xc179,0x0fb0,0x7fa6,
    0x3fd2,0x6b10,0xa2eb,0x45a8,
    0x3fb2,0x2755,0x5046,0x8d61,
    0x3f7a,0x7c93,0x2883,0x71d4,
    0x3f30,0x0ecf,0x118c,0x37ac,
    0x3ed0,0xe8f4,0xfcf7,0x4b5c,
    0x3e5e,0xf7f1,0x9edf,0xce99,
    0x3dd7,0xdc0d,0xdc57,0xc710,
    0x3d3c,0xe0fd,0x54c8,0xcd17,
    0x3c88,0x677f,0x1c0f,0x7fe8,
    0x3bb7,0x6b92,0xc6f8,0x96ab,
    0x3ac0,0x6494,0x9f50,0xaa8e,
    };
    static short B4[48] = {
      /*0x3ff0,0x0000,0x0000,0x0000,*/
    0x4006,0x91f2,0x05e7,0x2e77,
    0x3ffc,0x7bc9,0x2570,0x37d6,
    0x3fd8,0x2f54,0x9821,0x5c12,
    0x3fa0,0x5a4a,0xa7af,0x4bce,
    0x3f53,0x06a5,0x4a3a,0xfc75,
    0x3ef3,0x94a9,0xa71a,0x7eb2,
    0x3e81,0xb420,0xb4f6,0x00ce,
    0x3dfb,0x182b,0x1c44,0xba10,
    0x3d60,0x573c,0xe32c,0x40a9,
    0x3cab,0x92f6,0xc77b,0x077b,
    0x3bda,0x721c,0x6851,0xe1fe,
    0x3ae2,0x81e1,0x5f6f,0xbc13,
    };
    #endif
    
    
    #ifdef ANSIPROT
    extern double spence ( double );
    extern double polevl ( double, void *, int );
    extern double p1evl ( double, void *, int );
    extern double zetac ( double );
    extern double md_pow ( double, double );
    extern double md_powi ( double, int );
    extern double md_log ( double );
    extern double fac ( int i );
    extern double md_fabs (double);
    double polylog (int, double);
    #else
    extern double spence(), polevl(), p1evl(), zetac();
    extern double md_pow(), md_powi(), md_log();
    extern double fac(); /* factorial */
    extern double md_fabs();
    double polylog();
    #endif
    extern double MACHEP;
    
    
    double
    polylog (n, x)
         int n;
         double x;
    {
      double h, k, p, s, t, u, xc, z;
      int i, j;
    
    
    /*  This recurrence provides formulas for n < 2.
    
    
        d                 1
        --   Li (x)  =   ---  Li   (x)  .
        dx     n          x     n-1
    
    
    */
    
    
      if (n == -1)
        {
          p  = 1.0 - x;
          u = x / p;
          s = u * u + u;
          return s;
        }
    
    
      if (n == 0)
        {
          s = x / (1.0 - x);
          return s;
        }
    
    
      /* Not implemented for n < -1.
         Not defined for x > 1.  Use cpolylog if you need that.  */
      if (x > 1.0 || n < -1)
        {
          mtherr("polylog", DOMAIN);
          return 0.0;
        }
    
    
      if (n == 1)
        {
          s = -md_log (1.0 - x);
          return s;
        }
    
    
      /* Argument +1 */
      if (x == 1.0 && n > 1)
        {
          s = zetac ((double) n) + 1.0;
          return s;
        }
    
    
      /* Argument -1.
                            1-n
         Li (-z)  = - (1 - 2   ) Li (z)
           n                       n
       */
      if (x == -1.0 && n > 1)
        {
          /* Li_n(1) = zeta(n) */
          s = zetac ((double) n) + 1.0;
          s = s * (md_powi (2.0, 1 - n) - 1.0);
          return s;
        }
    
    
    /*  Inversion formula:
     *                                                   [n/2]   n-2r
     *                n                  1     n           -  md_log    (z)
     *  Li (-z) + (-1)  Li (-1/z)  =  - --- md_log (z)  +  2  >  ----------- Li  (-1)
     *    n               n              n!                -   (n - 2r)!    2r
     *                                                    r=1
     */
      if (x < -1.0 && n > 1)
        {
          double q, w;
          int r;
    
    
          w = md_log (-x);
          s = 0.0;
          for (r = 1; r <= n / 2; r++)
        {
          j = 2 * r;
          p = polylog (j, -1.0);
          j = n - j;
          if (j == 0)
            {
              s = s + p;
              break;
            }
          q = (double) j;
          q = md_pow (w, q) * p / fac (j);
          s = s + q;
        }
          s = 2.0 * s;
          q = polylog (n, 1.0 / x);
          if (n & 1)
        q = -q;
          s = s - q;
          s = s - md_pow (w, (double) n) / fac (n);
          return s;
        }
    
    
      if (n == 2)
        {
          if (x < 0.0 || x > 1.0)
        return (spence (1.0 - x));
        }
    
    
    
    
    
    
      /*  The power series converges slowly when x is near 1.  For n = 3, this
          identity helps:
    
    
          Li (-x/(1-x)) + Li (1-x) + Li (x)
            3               3          3
                         2                               2                 3
           = Li (1) + (pi /6) md_log(1-x) - (1/2) md_log(x) md_log (1-x) + (1/6) md_log (1-x)
               3
      */
    
    
      if (n == 3)
        {
          p = x * x * x;
          if (x > 0.8)
        {
          u = md_log(x);
          s = u * u * u / 6.0;
          xc = 1.0 - x;
          s = s - 0.5 * u * u * md_log(xc);
              s = s + PI * PI * u / 6.0;
              s = s - polylog (3, -xc/x);
          s = s - polylog (3, xc);
          s = s + zetac(3.0);
          s = s + 1.0;
          return s;
        }
          /* Power series  */
          t = p / 27.0;
          t = t + .125 * x * x;
          t = t + x;
    
    
          s = 0.0;
          k = 4.0;
          do
        {
          p = p * x;
          h = p / (k * k * k);
          s = s + h;
          k += 1.0;
        }
          while (md_fabs(h/s) > 1.1e-16);
          return (s + t);
        }
    
    
    if (n == 4)
      {
        if (x >= 0.875)
          {
        u = 1.0 - x;
        s = polevl(u, A4, 12) / p1evl(u, B4, 12);
        s =  s * u * u - 1.202056903159594285400 * u;
        s +=  1.0823232337111381915160;
        return s;
          }
        goto pseries;
      }
    
    
    
    
      if (x < 0.75)
        goto pseries;
    
    
    
    
    /*  This expansion in powers of md_log(x) is especially useful when
        x is near 1.
    
    
        See also the pari gp calculator.
    
    
                          inf                  j
                           -    z(n-j) (md_log(x))
        polylog(n,x)  =    >   -----------------
                           -           j!
                          j=0
    
    
          where
    
    
          z(j) = Riemann zeta function (j), j != 1
    
    
                                  n-1
                                   -
          z(1) =  -md_log(-md_log(x)) +  >  1/k
                                   -
                                  k=1
      */
    
    
      z = md_log(x);
      h = -md_log(-z);
      for (i = 1; i < n; i++)
        h = h + 1.0/i;
      p = 1.0;
      s = zetac((double)n) + 1.0;
      for (j=1; j<=n+1; j++)
      {
        p = p * z / j;
        if (j == n-1)
          s = s + h * p;
        else
          s = s + (zetac((double)(n-j)) + 1.0) * p;
      }
      j = n + 3;
      z = z * z;
      for(;;)
        {
          p = p * z / ((j-1)*j);
          h = (zetac((double)(n-j)) + 1.0);
          h = h * p;
          s = s + h;
          if (md_fabs(h/s) < MACHEP)
        break;
          j += 2;
        }
      return s;
    
    
    
    
    pseries:
    
    
      p = x * x * x;
      k = 3.0;
      s = 0.0;
      do
        {
          p = p * x;
          k += 1.0;
          h = p / md_powi(k, n);
          s = s + h;
        }
      while (md_fabs(h/s) > MACHEP);
      s += x * x * x / md_powi(3.0,n);
      s += x * x / md_powi(2.0,n);
      s += x;
      return s;
    }

  5. #5
    and the hat of wrongness Salem's Avatar
    Join Date
    Aug 2001
    Location
    The edge of the known universe
    Posts
    32,553
    Seems an awful lot like you just found some code, and you're still stuck wondering what to do with it all.
    http://cpansearch.perl.org/src/RKOBE...ibmd/polylog.c
    If you dance barefoot on the broken glass of undefined behaviour, you've got to expect the occasional cut.
    If at first you don't succeed, try writing your phone number on the exam paper.
    I support http://www.ukip.org/ as the first necessary step to a free Europe.

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Replies: 2
    Last Post: 12-10-2012, 11:25 PM
  2. Replies: 1
    Last Post: 03-03-2009, 03:47 PM
  3. Replies: 5
    Last Post: 08-16-2007, 11:43 PM
  4. Replies: 18
    Last Post: 11-13-2006, 12:11 PM
  5. making a program leave a msg for background program when it closes
    By superflygizmo in forum Windows Programming
    Replies: 2
    Last Post: 02-06-2006, 06:44 PM

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21