integrating discontinuous function

Hi there, I want to wirte a program which can integrating a defined function as follows

function f=r(x,p) if x<=0 f=0; else f=x^p; end

and use sybolic integration tool to integrate defined function

int(r(y-a),a,x));

the problem is matlab does not know whether to select f=0 or f=x^p.

please show me a way to solve the problem

Thanks

L.

Answer 1

[Andrew Newell]

Did you define your symbolic variables first?

syms a x y p

You'll need to give a numeric value for p or else int won't be able to evaluate it:

p=2;

You can redefine your function in terms of a heaviside function. Now you can integrate it:

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

fint =

    -(a - x)^3/3

If you want to get an answer for a given a and x, you can substitute:

subs(fint,{a,x},{sym(-1),sym(1)})

ans =

    8/3

double(ans)

ans =

    2.6667

EDIT: As Susan points out, this gives the wrong answer for a<x. This seems to occur because one of the integral limits, a, is also a parameter in the equation. This code gives the correct answer:

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

ans =

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

subs(fint,{x,a},{sym(2),sym(1)})

ans =

1/3

subs(fint,{x,a},{sym(1),sym(2)})

ans =

     0

[Lam]

Thanks Andrew. Yes I did define my symbolic variables. The problem is when integrating int(r(y-a,2),a,x)) with y is dummy variable. My final purpose is to plot the resulting integral when x varies from -1 to 1.

[Lam]

heaviside(x) gives you 0.5 when x==0. In my case, f=0 if x<=0.

[Andrew Newell]

That won't affect any integral. It's a basic fact of calculus that you can change any function at a single point (or even an infinite number of points, as long as they aren't "dense"), and have no effect on the integral.

[Susan]

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?

syms a x y p

p=2;

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

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

fint =

-(a - x)^3/3

ans =

-9

[Andrew Newell]

Thanks for pointing that out, Susan. I have corrected the answer.