Integrating Absolute Values

Doug Meade's picture

Here's an interesting problem from a colleague. Consider the following two improper integrals; the only difference is the absolute value in the integrand. I think this is another instance in which Maple mis-applies the FTOC.

restart;
with( IntegrationTools ):
#infolevel[all]:=3:
q := x*sin(1/x^2)-cos(1/x^2)/x;
                                            /1 \
                                         cos|--|
                                            | 2|
                                  /1 \      \x /
                             x sin|--| - -------
                                  | 2|      x   
                                  \x /          
I0 := Int( q, x=0..1 ):
I1 := Int( abs(q), x=0..1 ):

Maple gives values for both integrals.


value( I0 );
                                  1       
                                  - sin(1)
                                  2       
value( I1 );
                                  1       
                                  - sin(1)
                                  2       

It's a little odd that these two answers came out the same. It's easy to see that q does change signs infinitely often in [0,1]. A nice visualization of this can be obtained with:


plot([1/x,abs(q)],x=0..1,y=0..50);

We can dig a little deeper by looking at the values of the definite integrals on [A,1] and taking a limit as A->0.


int( q, x=A..1 );
                            1  2    /1 \   1       
                          - - A  sin|--| + - sin(1)
                            2       | 2|   2       
                                    \A /           
int(abs(q),x=A..1);
                        / 2    /1 \      /1 \\                   
                        |A  sin|--| - cos|--||                   
                        |      | 2|      | 2||                   
             1  2       |      \A /      \A /|    /1 \   1       
           - - A  signum|--------------------| sin|--| + - sin(1)
             2          \         A          /    | 2|   2       
                                                  \A /           

This starts to show where Maple is missing some important information. More is seen by looking at the antiderivatives of q and abs(q):


Q := int(q,x);
                                1  2    /1 \
                                - x  sin|--|
                                2       | 2|
                                        \x /
Q2 := int( abs(q), x );
                             / 2    /1 \      /1 \\        
                             |x  sin|--| - cos|--||        
                             |      | 2|      | 2||        
                  1  2       |      \x /      \x /|    /1 \
                  - x  signum|--------------------| sin|--|
                  2          \         x          /    | 2|
                                                       \x /

plot( [q, Q], x=0..1, view=[DEFAULT,-5..5] ):

plot( [abs(q), Q2], x=0..1, view=[DEFAULT,-1..5] );

plot( [Q,Q2], x=0..0.6, thickness=[3,1], discont=true );

Clearly, the two improper integrals at the start of this post should not have the same values.

Let's convert the problem to an improper integral on [1,infinity) by making a substitution x=1/u. When Maple makes this change of variable, it finds:


I2 := Change( I1, x=1/u );
                      /infinity       / 2\      / 2\  2   
                     |            -sin\u / + cos\u / u    
                   - |          - --------------------- du
                     |                      3             
                    /1                     u              

It's a little interesting that the absolute values have disappeared. (A plot quickly shows that this integrand is not always an increasing function.


q2 := GetIntegrand( Combine(I2) );
                                / 2\      / 2\  2
                            -sin\u / + cos\u / u 
                            ---------------------
                                      3          
                                     u           
plot( q2, u=1..50 );

What gets really interesting is the value of this integral:

 

value( I2 );
                                   1       
                                 - - sin(1)
                                   2       

Maple returns the same value for the definite integral with the correct (absolute value) in the integrand. Once again it is surprising that Maple reports the same values for these two integrals. But, the fact that these values are negative - even when the integrand is non-negative - is alarming.

It does appear that in some cases Maple is able to detect and correctly handle an infinite number of discontinuities in the antiderivative:


int(abs(cos(v)/v),v=1..infinity);
                /1   \
     -Ci(1) + Ci|- Pi|
                \2   /

          /infinity                                                   \
          | -----                                                     |
          |  \                                                        |
          |   )    /  /1           \     _k     /3           \     _k\|
        + |  /     |Ci|- Pi + Pi _k| (-1)   - Ci|- Pi + Pi _k| (-1)  ||
          | -----  \  \2           /            \2           /       /|
          \ _k = 0                                                    /
evalf( % );
                               Float(infinity)

 

 

Axel Vogt's picture

'Change'

only a minor comment:

  PDEtools[dchange](x=1/u ,I1);
  student[changevar](x=1/u ,I1,u);
seem not to ignore the absolute value.

check sign?

I wonder whether the code for definite integration includes a check on whether the integrand does not change sign within the interval of integration.

For  an integrand  that  is an 'abs'  call  and some simple variations it  is trivial.  But  for more complex  functions  it could be  difficult  (if not  impossible).  For this reason it might be  that there is no such check.

 

Axel Vogt's picture

Bug in 'simplify' or 'abs' ?

May be that comes through a bug in 'simplify' or 'abs'.

Define I3:=PDEtools[dchange](x=1/u ,I1) with integrand p:=op(1,%) =
abs(-1/u*sin(u^2)+cos(u^2)*u) / u^2, which is your I2, but with the
absolute sign is not forgotten.

Now apply 'simplify' to I3:

  simplify(I3);

                    infinity
                   /               2         2   2
                  |          -sin(u ) + cos(u ) u
                  |          --------------------- du
                  |                    3
                 /                    u
                   1

which Maple evaluates to -1/2*sin(1) and that coincides with the
given value for I0 ( u=1/x shows it to - I0).



This leads to find the following bug: 

Both simplify(p) and simplify(p,symbolic) behave correct. However
they do not if using assumptions:

  q;


                      |        2              |
                      |   sin(u )        2    |
                      | - ------- + cos(u ) u |
                      |      u                |
                      -------------------------
                                  2
                                 u

  simplify(q) assuming 0<u;
  # plot(%,u=1 .. sqrt(2*Pi));

                              2         2   2
                        -sin(u ) + cos(u ) u
                        ---------------------
                                  3
                                 u

To be sure that this is not caused through the 'assuming' command
one produces the same error through 'assume':

  assume(0<v): eval(p, u=v): simplify(%);

                              2         2   2
                        -sin(v ) + cos(v ) v
                        ---------------------
                                  3
                                 v

Strange that eval(p, u=sqrt(v)), eval(p, u=v^2) or eval(p, u=(v+1))
work as expected.

Since it is not unlikely that IntegrationTools[Change] uses 'simplify'
that would 'explain' it and would not mean that the integration itself
has an error here: may be 'Change' wants to use information about the
range of integration.


In fact the reason might be in the 'abs' function:

  r:=(sin(v^2)+cos(v^2)*v^2)/v^3; # we have 0 < v
  s:=abs(r);
                                2         2   2
                           sin(v ) + cos(v ) v
                      r := --------------------
                                     3
                                    v


                                2         2   2
                           sin(v ) + cos(v ) v
                      s := --------------------
                                     3
                                    v

Which is false, just plot it. While eval(r, v=u); abs(%); works ok ...
JacquesC's picture

The bug is in 'is', possibly `is/RealRange`

with r defined as just sin(v^2)+cos(v^2)*v^2 (the numerator) under the assumption v>0 , we get

is(0 < r), is(0 <= r);
                             false, true

where the second one is clearly a bug.  In fact, as far as I can tell, the bug is in `is/RealRange` the line

`property/PropInclusion`(objprop,newprop)

has things backwards, ie it should really be

 

`property/PropInclusion`(newprop,objprop)

If I am right, I wonder how many bugs in is (and thus signum, etc) are due to this?

thanks

We can look into this. Thanks, Jacques.

Dave Linder
Mathematical Software, Maplesoft

work to be done

Currently,

> assume(n>0):
> is(2^n>n);
                                    true

> assume(r,RealRange(Open(-1),Open(1)));
> limit(r^n,n = infinity);
                                            n
                                lim       r~
                           n -> infinity

And, after the simple flip of how `is/RealRange` calls `property/PropInclusion` in its last line,

> assume(n>0):
> is(2^n>n);
                                     false

> assume(r,RealRange(Open(-1),Open(1)));
> limit(r^n,n = infinity);
                                       0

There's always work to be done.

Dave Linder
Mathematical Software, Maplesoft

JacquesC's picture

refined guess

I have a refinement on what I think might be going on here.  It is possible that what is wrong is not `is/RealRange` itself but rather how it is being called.  When tracing things, what was passed to this routine and how it was processed looked a tad odd.  Maybe in some case it is being called with one convention, while in another it is using an opposite convention?

I have not looked at this in detail, this is just a guess.  It would take me a couple of hours to look at that carefully, a couple of hours I should really spend doing something else... [even though tracing this bug might just be more fun],

where is the bug?

See this also:

is(s>0);
                                 true

PD: Jacques posted while I was editing.

Axel Vogt's picture

thx for nailing it down

thx for nailing it down, that was worth the effort

Axel Vogt's picture

FWIW: numerical

just want to point out, that these days there was a quite similar task at the Maple NG, where the
poster was looking for the numerical value (but never said why)

groups.google.de/group/comp.soft-sys.math.maple/browse_thread/thread/b3710ca0dc842857/#

Comment viewing options

Select your preferred way to display the comments and click "Save settings" to activate your changes.
}