Re: Integrating x^b*Log ^m gives wrong result?

*To*: mathgroup at smc.vnet.net*Subject*: [mg85449] Re: [mg85434] Integrating x^b*Log ^m gives wrong result?*From*: danl at wolfram.com*Date*: Mon, 11 Feb 2008 06:05:25 -0500 (EST)*References*: <200802101019.FAA17904@smc.vnet.net>

> Dear all, > > I'm running into the following problems with symbolic vs. numerical > integration of the function x^(-3.5)*Log[x]^m: > > In[564]:= > ClearAll["Global`*"]; > f1[m_]:=N[Integrate[x^(-3.5)*Log[x]^m,{x,5,10}]]; > f2[m_]:=NIntegrate[x^(-3.5)*Log[x]^m,{x,5,10}]; > Map[f1,{5,10,25,40}] > Map[f2,{5,10,25,40}] > > Out[567]= {0.145434,4.62609,401145.,-9.30763*10^23} > Out[568]= {0.145434,4.62609,403156.,6.33616*10^10} > > Of course the symbolic integration is wrong here since it shouldn't > yield a negative number. If the recursive formula resulting from > partial integration is used, things seem to go wrong as well: > > In[572]:= > f[m_]:=(-1/2.5)*(10^(-2.5)*Log[10]^m-5^(-2.5)*Log[5]^m)+(m/ > 2.5)*f[m-1]; > f[0]=(-1/2.5)*(10^(-2.5)-5^(-2.5)); > Map[f,{5,10,25,40}] > > Out[574]= {0.145434,4.62609,403156.,-2.54037*10^16} > > So the result for m=25 still coincides with the one from NIntegrate, > while Integrate already gives something different; for m=40 the result > is different from both NIntegrate and Integrate (and wrong as it is > negative). If one changes the negative power of x to a positive one, > things seem ok btw. > > Any clues what might be going on here? > > Thanks in advance, Kees > > In[533]:= $Version > Out[533]= 6.0 for Microsoft Windows (32-bit) (April 27, 2007) > The issue is one of massive cancellation error. Here I use "error" in the sense of numerical analysis, and not to indicate a bug. If you take the indefinite integral for, say, m=40, you will see many terms with large coefficients, and both positive and negative contributions. When evaluated at the interval endpoints these will tend to cancel in the significant digits. At m=40, all meaningful digits have been lost. How to avoid this? You will need to do two things. (1) Use exact arithmetic throughout the integration. So change to f1[m_] := Integrate[x^(-7/2)*Log[x]^m, {x, 5, 10}] (2) When evaluating numerically, use bignum rather than default machine arithmetic. Say N[f1[40],20]. The internals of N will either figure out the precision needed to attain 20 digits in the result, or else warn you of failure to accomplish that. As a general remark, mixing approximate arithmetic with symbolic computation is a dicey business. Numeric evaluation of exact results is (almost always) a reasonable thing to do. But first obtain the exact result, or else use dedicated numeric methods (such as NIntegrate) from the start. Daniel Lichtblau Wolfram Research

**References**:**Integrating x^b*Log[x]^m gives wrong result?***From:*KvS <keesvanschaik@gmail.com>