Integrating piecewise function

I'm following up on this q&a:

http://www.mathworks.com/matlabcentral/answers/3653-integrating-discontinuous-function

I typed in the following code and I don't understand the result:

syms a x y p
p=2;
fint = int(heaviside(y-a)*(y-a)^p,y,a,x)
subs(fint,{a,x},{sym(0),sym(-3)})

fint =
-(a - x)^3/3

ans =
-9

I guess I'm being dense here, but it seems that the heaviside function has had no effect. Shouldn't answers for x<a be zero? Why isn't fint a piecewise function?

Thanks,

Susan

[Andrew Newell]

I have corrected the original answer.

Answer 1

[Walter Roberson]

It would be easier to read (and might even solve the problem) if you were to explicitly specify the variable of integration in the int() call instead of just specifying the bounds and having it try to deduce the variable.

[Susan]

Ok, I edited to specify y as variable of integration. There is no change in result.

[Susan]

Interestingly, the following change gets me the right answer (though it's not as general as the original case):

p=2;

fint = int(heaviside(y-a)*(y-a)^p,y,0,x)

subs(fint,{a,x},{sym(0),sym(-3)})

fint =

(a^3*heaviside(-a))/3 - (heaviside(x - a)*(a - x)^3)/3

ans =

0

So there is some problem with simultaneously using a in the heaviside argument and as a limit of integration. The original result still seems wrong to me though.

[Walter Roberson]

Just to check: are you using the Maple based symbolic toolbox (older versions) or the MuPad based one (newer versions) ?

When I try this out in Maple directly, the answer I get after substitution is

limit(-Heaviside(y)*y^(p+1)/(p+1)+3^(p+1)/(p+1), y = 0, right)

which simplifies to

-(limit(Heaviside(y)*y^(p+1)-3^(p+1), y = 0, right))/(p+1)

The piecewise() equivalent of this is (in Maple notation)

piecewise(y < 0, 3^(p+1)/(p+1), y = 0, -(limit(-3^(p+1)+y^(p+1)*undefined, y = 0, right))/(p+1), 0 < y, -(limit(-3^(p+1)+y^(p+1), y = 0, right))/(p+1))

Maple's "undefined" is pretty much equivalent to NaN -- that is, Heaviside is not defined when its argument is 0.

Answer 2

[Susan]

So it looks to me that, if I'm careful, I can get the right answer:

%%what does fint(x) look like for fixed a? 
% define a separate lower limit of integration
syms a x y p x0
p=2; % power
a=1; % function is zero-valued below a
x0=0;% start integrating here
fint = int(heaviside(y-a)*(y-a)^p,y,x0,x)

% ezplot(fint)  % doesn't work for piecewise syms?
xx = linspace(-3,6,100);
yy = subs(fint,x,xx);
plot(xx,yy);

% great, exactly what I think is right
fint =
piecewise([x < 1, 0], [1 <= x, (heaviside(x - 1)*(x - 1)^3)/3])

But for some lower limits of integration, the need for a piecewise solution is missed entirely:

%%what if x0 ==2?

x0=2;% start integrating here
fint = int(heaviside(y-a)*(y-a)^p,y,x0,x)

% ezplot(fint)  % doesn't work for piecewise syms?
xx = linspace(-3,6,100);
yy = subs(fint,x,xx);
plot(xx,yy);

% eww.  Wrong answer!
fint =
((x - 2)*(x^2 - x + 1))/3

[Walter Roberson]

If you mentally add the assumption that x >= x0, then when you use an x0 that is greater than your a, then x < a is going to be false so you fall over to just the heaviside portion.

[Susan]

So this assumption is obvious? I guess the lesson is that mupad will make assumptions and the user needs to be wary of them.

Answer 3

[Christopher Creutzig]

When integrating from a to x, int makes the implicit assumption that a ≤ y ≤ x and thus a ≤ x, unless that is obviously wrong, in which case the interval borders are swapped. This is a side-effect of restricting the variable of integration to the interval given.

[Susan]

Same comment as above to Walter: I guess the lesson is that mupad will make assumptions and the user needs to be wary of them.