Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
add comment

3 Answers

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}]]

plot

share|improve this answer
add comment

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}}]]

enter image description here

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).

share|improve this answer
add comment

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.

share|improve this answer
    
Thank you for your suggestions. I will rewrite the definitions. –  Wayan Apr 3 '12 at 1:36
    
@wayan I'm not so sure you should replace Boole with Piecewise. As the comments under this answer show, MMA does smart things when Boole and NIntegrateare combined. Maybe it does here too. –  Sjoerd C. de Vries Apr 3 '12 at 10:53
    
@SjoerdC.deVries Yes, I try indicator[a_,b_][x_] as recommended by @R.M but the results are different. Still working on it. –  Wayan Apr 3 '12 at 12:04
    
@SjoerdC.deVries I'm not fully convinced by Heike's example there. It probably has to do with the way the OP constructed Piecewise vs. Boole and wrong constructions can lead to different results. For example, replace Boole[x + y < .5] f[x, y] in Heike's example with g[x, y], where g[x_, y_] := Piecewise[{{f[x, y], x + y < 0.5}, {0, True}}], and you'll see the same set of points being returned with the same time taken. –  rm -rf Apr 3 '12 at 15:15
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.