program

[Ubuntu2014]

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.

[Salem]

Well the first thing you should do is post your current attempt.

[grumpy]

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.

[Ubuntu2014]

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 :

/*                            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;
}

[Salem]

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