Integrating with multiple indicator functions

Question

I am trying to calculate an integral involving multiple indicator functions, such as:

$$ h(u,v,w) = -\int_0^1 J^{\prime\prime}(s) (I_{(0,s]}(u) - s)(I_{(0,s]}(v) - s)(I_{(0,s]}(w) - s)\, dF^{-1}(s)$$

where $0 < u, v, w < 1$ and

$$\begin{align} J(s) &= 6s(1-s)\\ F(s) &= 1/[1 + \exp[-s]]\\ F^{-1}(s) &= -\log[-1 + 1/s]\\ \end{align} $$

Here is my code:

Js = 6 s (1 - s) (*weight function*)
FDJs = D[Js, s] (*first derivative of Js*)
SDJs = D[FDJs, s] (*second derivative of Js*)
Fs = 1/(1 + Exp[-(s)]) (*Fs*)
Finv = -Log[-1 + 1/s] (*F inverse*)
DFinv = D[Finv, s] (*derivative of F inverse*)
h = -Integrate[
    SDJs*(Boole[u <= s] - s)*(Boole[v <= s] - s)*(Boole[w <= s] - 
    s) DFinv, {s, 0, 1}, 
    Assumptions -> {0 < u < 1 , 0 < v < 1, 0 < w < 1}]

I get some output, but I am not sure if the syntax is correct especially for the indicator function.

Answer 1

A better construct for your indicator function would be using Piecewise. So you'd define an indicator function $I_{(a,b]}$ as:

indicator[a_, b_][x_] := Piecewise[{{1, a < x <= b}, {0, True}}]

However, there are other issues in your code that could be improved.

• First, avoid using capital letters for your variables/functions. This is because Mathematica uses that notation for its own functions and invariably, you'll (at some point) be bitten by it.
• Define your functions by appropriately scoping the variable. For example you defined Js = 6 s (1 - s). However, if s had a pre-existing value (you never know... sometimes it could be in a different notebook but shared context), then you're assigning to Js the result based on that value of s.
• Use Set (=) instead of SetDelayed (:=) when you use D. Otherwise, you'll end up calculating the derivative at each integration point. See this answer for an example of why this matters.

My suggested rewrite of the definitions part of your problem would be something like this:

j[s_] := 6 s (1 - s);      (*weight function*)
j1[s_] = D[j[s], s] ;      (*first derivative of Js*)
j2[s_] = D[j[s], {s, 2}];  (*second derivative of Js*)
f[s_] := 1/(1 + Exp[-(s)]) (*Fs*)
finv[s_] := -Log[-1 + 1/s] (*F inverse*)
finv1[s_] = D[finv[s], s]; (*derivative of F inverse*)

You can try working with this and the indicator function above to see if it makes any difference or if you're getting the answer you expect.

Answer 2

Assuming definitions of the functions as in R.M.'s answer we can observe that :

{ j1[s], j2[s], finv1[s] } // Simplify

{6 - 12 s, -12, 1/(s - s^2)}

To determine the integrand we need the functions : $ I_{(0,s]}(u) - s \quad$ and $ (I_{(0,s]}(u) - s)(I_{(0,s]}(v) - s)(I_{(0,s]}(w) - s)$, so we define them this way :

indicat[a_, b_, x_] := UnitStep[x - a] UnitStep[b - x] - b

indi[u_, v_, w_, s_] := indicat[0, s, u] indicat[0, s, v] indicat[0, s, w]

A more reasonable idea is to use NIntegrate instead of Integrate and to make use of j2[s] and finv1[s] to define the function h :

h[u_, v_, w_] :=  12 NIntegrate[1/(s - s^2) indi[u, v, w, s], {s, 0, 1}]

It appears that Integrate involves some drawbacks, (see e.g. this answer) and it would take much more time. Now we can make 3D-plots of h for various arguments, e.g. w :

Grid[ Table[ Plot3D[ h[u, v, w], {u, 0, 1}, {v, 0, 1}], {w, {0.2, 0.4, 0.6, 0.8}}]]

It takes rather a long time to make these plots working with definitions of indicators by Piecewise as well as by UnitStep(approximately the same).

Answer 3

One can calculate the integral observing that h[u,v,w] is symmetric in it arguments and that one can split the integration interval in [0,u]\[Union][u,v]\[Union][v,w]\[Union][w,1]. In each interval the indicator is either 0 or 1. The result is :

12 Simplify[Integrate[1/(s - s^2) (UnitStep[u] UnitStep[s - u] - s) (UnitStep[v] UnitStep[s - v] - s) (UnitStep[w] UnitStep[s - w] - s), {s, 0, u}, Assumptions -> {0 < u < 1, 0 < v < 1, 0 < w < 1, u < v < w}] + Integrate[1/(s - s^2) (UnitStep[u] UnitStep[s - u] - s) (UnitStep[v] UnitStep[s - v] - s) (UnitStep[w] UnitStep[s - w] - s), {s, u, v}, Assumptions -> {0 < u < 1, 0 < v < 1, 0 < w < 1, u < v < w}] + Integrate[1/(s - s^2) (UnitStep[u] UnitStep[s - u] - s) (UnitStep[v] UnitStep[s - v] - s) (UnitStep[w] UnitStep[s - w] - s), {s, v, w}, Assumptions -> {0 < u < 1, 0 < v < 1, 0 < w < 1, u < v < w}] + Integrate[1/(s - s^2) (UnitStep[u] UnitStep[s - u] - s) (UnitStep[v] UnitStep[s - v] - s) (UnitStep[w] UnitStep[s - w] - s), {s, w, 1}, Assumptions -> {0 < u < 1, 0 < v < 1, 0 < w < 1, u < v < w}]]

(* 12 (-(3/2) + u + v + w + Log[1 - u] - Log[w]) *)

hAux[u_?NumericQ,v_?NumericQ,w_?NumericQ]:=12 (-(3/2)+u+v+w+Log[1-u]-Log[w])
h[u_?NumericQ, v_?NumericQ, w_?NumericQ] := hAux[Sequence @@ Sort[{u, v, w}]]