Integrating piecewise function

5 次查看（过去 30 天）
Susan
Susan 2011-5-25
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

请先登录，再回答此问题。

回答（3 个）

Walter Roberson
Walter Roberson 2011-5-25
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.
  3 个评论
显示 1更早的评论
Susan
Susan 2011-5-25
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
Walter Roberson 2011-5-25
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.

请先登录，再进行评论。

Susan
Susan 2011-5-25
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
  2 个评论
Walter Roberson
Walter Roberson 2011-5-25
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
Susan 2011-6-2
So this assumption is obvious? I guess the lesson is that mupad will make assumptions and the user needs to be wary of them.

请先登录，再进行评论。

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.
  1 个评论
Susan
Susan 2011-6-2
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Numeric Solvers 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by