       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:=
> 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= {0.145434,4.62609,401145.,-9.30763*10^23}
> Out= {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:=
> f[m_]:=(-1/2.5)*(10^(-2.5)*Log^m-5^(-2.5)*Log^m)+(m/
> 2.5)*f[m-1];
> f=(-1/2.5)*(10^(-2.5)-5^(-2.5));
> Map[f,{5,10,25,40}]
>
> Out= {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?
>
>
> In:= \$Version
> Out= 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,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

```

• Prev by Date: Re: Import syntax question
• Next by Date: Re: Where to put ConvexHull3D.m
• Previous by thread: Integrating x^b*Log[x]^m gives wrong result?
• Next by thread: Re: Integrating x^b*Log[x]^m gives wrong result?