Integral of piecewise functions

[Paul Butler]

Currently, taking the integral of a piecewise function in Sage gives
you the definite integral. I've proposed on trac that the integral of
piecewise functions be indefinite by default. This would be consistent
with how integration works on other functions in Sage, as well as
piecewise functions in Maple and Mathematica.

The main concern is whether the integral of a piecewise function is
even well-defined. It seems to me that at least for continuous
piecewise functions, the indefinite integral is well-defined. The
anti-derivative is well defined, and by the fundamental theorem of
calculus, the indefinite integral of a continuous function is the
anti-derivative. As for discontinuous piecewise functions, I'm finding
it difficult to convince myself either way.

The trac ticket is 4721 ( http://trac.sagemath.org/sage_trac/ticket/4721 )

-- Paul

[David Joyner]

First, the indefinite integral is the anti-derivative, by definition.
What the FTC says is that although the indefinite integral "evalated
at b" is not
well-defined, and the same antiderivative "evalated at a" is also not
well-defined,
their difference *is* well-defined. Moreover, this difference agrees with the
definite integral defined by the Riemann sum between a and b.

Second, I think (but I am not sure), what you want when you say
indefinite intergal is not the indefinite integral but is a function who
derivavtive is the original function, defined as follows:

if the orginial function f(x) is

f1(x), a1<x<=a2,
f2(x), a2<x<=a3,
...
fn(x), an<x<=a{n+1}

(and 0 outside (a1,a{n+1}) then I guess you want to define the integral,
call it F, by


int_{a1}^x f1(t) dt, a1<x<=a2,
int_{a2}^x f2(t) dt, a2<x<=a3,
...
int_{an}^x fn(t) dt, an<x<=a{n+1}

Is this correct? This is not the antiderivative
but it does have the property that F'(x)=f(x).

The antiderivative is only well-defined up to an additive constant.
IMHO, the piecewise defined function of antiderivavtives

int f1(x) dx +C1 , a1<x<=a2,
int f2(x) dx +C2, a2<x<=a3,
...
int fn(x) dx +Cn, an<x<=a{n+1}

does not make sense.

[Tim Lahey]

This crops up regularly in solid mechanics, in particular in the bending
of beams. One can represent the loading of a beam as derivatives and
integrals of Heaviside functions and then integrate the expression to
get the shear force, bending moment, and deflection of the beam. It's
also one of the parts of Maple that's problematic. I wrote some code
for my students a few years back to help them, but unfortunately, there
are problems with some integrals of Heaviside functions. These can be
discontinuous piecewise functions, at least the shear and moment graphs.

Cheers,

Tim.

---
Tim Lahey
PhD Candidate, Systems Design Engineering
University of Waterloo
http://www.linkedin.com/in/timlahey

[Paul Butler]

On Sat, Dec 6, 2008 at 12:37 PM, David Joyner <wdjo...@gmail.com> wrote:

On Sat, Dec 6, 2008 at 11:39 AM, Paul Butler <pau...@gmail.com> wrote:
>
> Currently, taking the integral of a piecewise function in Sage gives
> you the definite integral. I've proposed on trac that the integral of
> piecewise functions be indefinite by default. This would be consistent
> with how integration works on other functions in Sage, as well as
> piecewise functions in Maple and Mathematica.
>
> The main concern is whether the integral of a piecewise function is
> even well-defined. It seems to me that at least for continuous
> piecewise functions, the indefinite integral is well-defined. The
> anti-derivative is well defined, and by the fundamental theorem of
> calculus, the indefinite integral of a continuous function is the
> anti-derivative. As for discontinuous piecewise functions, I'm finding
> it difficult to convince myself either way.

First, the indefinite integral is the anti-derivative, by definition.
What the FTC says is that although the indefinite integral "evalated
at b" is not
well-defined, and the same antiderivative "evalated at a"  is also not
well-defined,
their difference *is* well-defined. Moreover, this difference agrees with the
definite integral defined by the Riemann sum between a and b.

Second, I think (but I am not sure), what you want when you say
indefinite intergal is not the indefinite integral but is a function who
derivavtive is the original function, defined as follows:

Actually, a function F which is simply the antiderivative of f is no use to me, I need a function with the property that F(b) - F(a) is the Riemann sum between a and b. (They are only guaranteed by the FTC to be the same function if f is continuous.) Maybe there is a better term for this?

if the orginial function f(x) is

f1(x), a1<x<=a2,
f2(x), a2<x<=a3,
...
fn(x), an<x<=a{n+1}

(and 0 outside (a1,a{n+1}) then I guess you want to define the integral,
call it F, by


int_{a1}^x f1(t) dt, a1<x<=a2,
int_{a2}^x f2(t) dt, a2<x<=a3,
...
int_{an}^x fn(t) dt, an<x<=a{n+1}

Is this correct? This is not the antiderivative
but it does have the property that F'(x)=f(x).

 
I think you and I are defining antiderivative differently. I'm using the definition that F is an antiderivative of f if F'(x) = f(x) for all x in the domain of f(x). (Also stated here: http://planetmath.org/encyclopedia/Antiderivative.html .)

Either way, the property F'(x) = f(x) is not necessarily true for piecewise antiderivatives defined that way. Consider this function.

f(x) = x, 0 <= x <= 1
f(x) = 1, 1 < x

If we use the definition you gave to find F = integral(f), F'(1) is undefined so it is not true that F'(x) == f(x) for all x.

Instead, we use the definition that F=

integrate(f1, t, a1, x), a1 < x <= a2
integrate(f2, t, a2, x) + integrate(f1, t, a1, a2), a2 < x <= a3
integrate(f3, t, a3, x) + integrate(f2, t, a2, a3) + integrate(f1, t, a1, a2), a3 < x <= a4
...
integrate(fn, t, an, x) + integrate(f[n-1], t, a[n-1], an) + ... + integrate(f1, t, a1, a2), an < x

(We also need a special case for when a1 = -infinity, which I didn't show.)

With this definition, F(b) - F(a) can be used to find the Riemann sum between a and b. Also, F'(x) = f(x) seems to hold, except at points where f(x) goes from defined to undefined or vice-versa.

The antiderivative is only well-defined up to an additive constant.
IMHO, the piecewise defined function of antiderivavtives

int  f1(x) dx +C1 , a1<x<=a2,
int  f2(x) dx +C2, a2<x<=a3,
...
int  fn(x) dx +Cn, an<x<=a{n+1}

does not make sense.

I agree that it doesn't make sense where C1 .. Cn are arbitrary constants.

-- Paul

[Ronan Paixão]

If I had the expertise to implement it, I would do the following:

The integration would return another Piecewise function in which the
first interval is integrated normally, and the next ones would integrate
the function in that interval and add the definite integral of the
previous intervals. I think that makes some sense.

Ronan Paixão

[David Joyner]

On Sat, Dec 6, 2008 at 3:56 PM, Ronan Paixão <ronan...@yahoo.com.br> wrote:
>
> If I had the expertise to implement it, I would do the following:
>
> The integration would return another Piecewise function in which the
> first interval is integrated normally, and the next ones would integrate
> the function in that interval and add the definite integral of the
> previous intervals. I think that makes some sense.

This makes sense for compactly supported functions.
How do you do that for something like f(x) = max(1,floor(x))?

[David Joyner]

On Sat, Dec 6, 2008 at 3:48 PM, Paul Butler <pau...@gmail.com> wrote:
>

>
> Either way, the property F'(x) = f(x) is not necessarily true for piecewise
> antiderivatives defined that way. Consider this function.
>
> f(x) = x, 0 <= x <= 1
> f(x) = 1, 1 < x
>
> If we use the definition you gave to find F = integral(f), F'(1) is
> undefined so it is not true that F'(x) == f(x) for all x.
>
> Instead, we use the definition that F=
>
> integrate(f1, t, a1, x), a1 < x <= a2
> integrate(f2, t, a2, x) + integrate(f1, t, a1, a2), a2 < x <= a3
> integrate(f3, t, a3, x) + integrate(f2, t, a2, a3) + integrate(f1, t, a1,
> a2), a3 < x <= a4
> ...
> integrate(fn, t, an, x) + integrate(f[n-1], t, a[n-1], an) + ... +
> integrate(f1, t, a1, a2), an < x
>
> (We also need a special case for when a1 = -infinity, which I didn't show.)

Okay, this helps me understand what you mean.

Still, the case a1=-ifinty is precisely the special case which I don't
understand.

For example, take a function such as f(x) = max(1,floor(x)), x real.
How do you define an antiderivative F(x) so that
F(b)-F(a) = area under the y=f(x) for a<x<b?
(And mayeb you can do it for that special function,
and let us ignore points of discontinuity for the sake of
discussion.) In other words, I am asking for the algorithmic procedure
you would use to create an "area function" of a piecewise-defined
function on the reals.

[Paul Butler]

When a1 = -infinity, I would make F1 = integrate(f, x, a2, x) instead of integrate(f, x, a1, x). Then I would not calculate the definite integral of the first interval, which would align my constants so that F(a2) = 0. When I get a chance, I'll add this to my code.

Functions like floor with an infinite number of pieces are beyond the scope of the Piecewise class, because piecewise functions in Sage can only (at the moment) have finitely many pieces. I can't give an algorithm for integration of all piecewise functions with multiple pieces off the top of my head, but I'll give it some thought.

-- Paul

[David Joyner]

I see what you want now. If you can make this work (ie, satisfy the FTC),
then it does make sense to use your code as the default.

Thanks for explaining your construction and for working on the code!

[Ronan Paixão]

> This makes sense for compactly supported functions.
> How do you do that for something like f(x) = max(1,floor(x))?
>

sage: f(x)=max(1,floor(x))
sage: f(10)
1
sage: type(f)
<class 'sage.calculus.calculus.CallableSymbolicExpression'>

huh, there are some problems with this function, first because it always
return 1 and second because it's not built with Piecewise, which is the
point of the discussion.

Ronan

[David Joyner]

Paul:
It was not at all clear to me when you first started working on the
antiderivatives of piecewise defined functions that it was possible or
made sense.
Now I have changed my mind to an extent. If you think of the
antiderivative as the
(multi-valued) inverse of the differentiation map, as I usually do, then I'm
not sure what to do. But if you just want a function whose derivavtive is the
original and which satisfies the FTC, then I think I have an algorithm which
might work.

Let f(x) be the piecewise defined function which is fn(x) for an<x<a(n+1),
for n in ZZ (and let us ignore endpots for now). Here a_n is an
infinite sequence
a_{-infty}=-infty and a_{infty}=+infty and each f_n(x) is Riemann integrable
with antiderivative Fn(x)+Cn. The question is can we piece the Cns's
together so that
F(x) = Fn(x)+Cn for an<x<a(n+1) satisfies the FTC? These constants must satisfy

Cn = F_{n-1}(an-)-Fn(an+)+C_{n-1},

If you fix C0 = C to be an arbitrary constant (you could take this to
be 0 for example),
this is a 1-step recursionrelation which can be solved.

I think this should do the job.

Thanks for raising this issue!

++++++++++++++++++++++++++++++++++++

[Paul Butler]

David, thanks.

What do you mean by an+ and an- in the notation?

From what I can tell, that does look like the same algorithm I plan to
use. I have some working code now, I just want to review it and make
sure I didn't miss any edge cases before re-posting the patch. The
code posted as a patch now uses the same algorithm, but doesn't
properly handle unbounded functions.

-- Paul

[David Joyner]

On Tue, Dec 9, 2008 at 9:30 PM, Paul Butler <pau...@gmail.com> wrote:
>
> David, thanks.
>
> What do you mean by an+ and an- in the notation?

In general, by c+ I mean the limit as you approach c from the right,
and c- is the limit as you approach c from the left.

>
> From what I can tell, that does look like the same algorithm I plan to
> use. I have some working code now, I just want to review it and make
> sure I didn't miss any edge cases before re-posting the patch. The
> code posted as a patch now uses the same algorithm, but doesn't
> properly handle unbounded functions.

Great. I discussed the algorithm with a few colleagues and I am pretty sure it
will work even if (when?) piecewise functions with infinitely many parts are
implemented in Sage. It's cool that once you fix one of the constants
of integration
they are all determined, thanks to the FTC, which is consistent with how the
inverse of the derivative map should behave.

Looking forward to your patch.