Integrating function with abs or max



On 1/22/07, Ray Tice <trayracing at yahoo.com> wrote:
> How do I reformulate functions that use abs or max so that maxima can
> integrate them?
> For example, integrate(abs(x), x, -1, 1)  won't evaluate, but both of these
> will:
>   integrate(abs(x), x, -1, 0)
>   integrate(abs(x), x,  0, 1)

Here is some code that handles the simpler cases.  It does *not*
handle nested max/min/abs, so your example of max(1-abs(x),0) will not
work.  It also fails silently when solve gives solutions like abs(x)=1
or x=abs(x).  I depends on the comparison package, which is rather
simple-minded.  But it does work for many cases (showing trace):

integrate_disc(abs(x)-abs(x+1),x,-2,2);
    integrate [abs(x) - abs(x + 1), x, - 2, - 1 => 1
    integrate [abs(x) - abs(x + 1), x, - 1, 0 => 0
    integrate [abs(x) - abs(x + 1), x, 0, 2 => - 2
=>    - 1

integrate_disc(max(x,x^2,x^3),x,-2,2);
  integrate max(x^3,x^2,x),x,-2,-1 => 7/3
  integrate max(x^3,x^2,x),x,-1,0 => 1/3
  integrate max(x^3,x^2,x),x,0,1 => 1/2
  integrate max(x^3,x^2,x),x,1,2 => 15/4
=>  83/12

integrate_disc(abs(x)/(x^2+max(1,x+1/2)),x,-1,2)
  integrate abs(x)/(max(x+1/2,1)+x^2),x,-1,0 => log(2)/2
  integrate abs(x)/(max(x+1/2,1)+x^2),x,0,1/2 => log(5/4)/2
  integrate abs(x)/(max(x+1/2,1)+x^2),x,1/2,2 =>
log(13)/2-atan(5)-log(5/2)/2+atan(2)
 => log(13)/2-atan(5)-log(5/2)/2+log(2)/2+atan(2)+log(5/4)/2

Comments, corrections, bug reports, and of course improvements are welcome.

                 -s

---------------------------------------------------
/* Find the breakpoints of an expression defined with abs or max */

inop(ex)   := inpart(ex,0)$
inargs(ex) := block([inflag:true], args(ex) )$

/* doesn't currently handle nested min/max/abs */
find_breaks(ex) :=
      if atom(ex)
	then {}
      elseif member(inop(ex),'[abs,signum]) then {inpart(ex,1)}
      elseif member(inop(ex),'[min,max]) then pairwise_differences(inargs(ex))
      elseif member(inop(ex),'[floor,ceiling]) then error("Can't
handle floor or ceiling")
      else xreduce('union, setify(map('find_breaks, inargs(ex))))$

find_zeroes(ex,var) :=
 block([r,s:[]],
   r:listify(xreduce(union,map(lambda([q],setify(map(rhs,solve(q,var)))),find_breaks(ex)),{})),
   while r#[] do
     /* ignore complex solutions */
     ( if asksign(imagpart(r[1]))=zero then s:cons(realpart(r[1]),s),
       r:rest(r) ),
     s)$

pairwise_differences(l) :=
   block( [ res: {} ],
	    while l # [] do
	       ( res: union(setify(makelist(if orderlessp(first(l)-i,i-first(l))
					       then i-first(l)
					       else first(l)-i,
					     i,rest(l))),
			    res),
		 l: rest(l) ),
	    res )$

sort_ask(l):=
  block([remove:{},s],
     sort(listify(l),
	  lambda([a,b],
		 s:asksign(a-b),
		 /* Just ignore this case
		 if s=zero then
		      remove: adjoin(if orderlessp(a,b) then [a,b] else [b,a],
				     remove),  */
		 if s=pos then false else true)))$

filter(f,l):=
  block([r:[]],
	while l#[] do
	   ( if not apply(f,[first(l)]) then r:cons(first(l),r), l:rest(l) ),
	reverse(r))$

integrate_disc(ex,x,low,high):=
  block([z,r:0],
    z: append([low],
	      sort_ask(filter(lambda([q],is(maybe(q<=low or q>=high)=true) ),
		              find_zeroes(ex,x))),
	      [high]),
    while length(z)>=2 do
      ( r:r+integrate(ex,x,z[1],z[2]),
	z: rest(z) ),
    r)$