Message Boards Message Boards

0
|
8116 Views
|
6 Replies
|
5 Total Likes
View groups...
Share
Share this post:

Integral evaluation fails

Welling Oei, UMC Utrecht
Posted 11 years ago

Hi All, I was trying to derive a formula by evaluating an integral. I defined all the restrictions for that setting and I wanted to have the solution in terms of a simplified formula. But Mathematica did not give the solution like last time (with some modification). I was wondering if there is anyone who can help me solve this problem.

Thank you in advance for your kind attention and support.

Here is the code for the equation A with Cp and infectable functions defined in it. 

Cp[te_, Di_, Dv_, D0_] :=
 Which[
  Di + D0 <= Dv + te || Dv + te <= 0, 0,
  D0 < Dv + te && D0 + te < 0 && Di > Dv + te, Inf/(Inf + Ip),
  D0 < Dv + te && D0 + te >= 0 && Di > Dv + te, (D0 Inf)/(
  D0 Inf - Ip te),
  Di <= Dv + 
     te && ((te == 0 && D0 >= Dv) || (D0 >= Dv + te && te >= 0)), Di/
  Dv,
  D0 >= Dv + te && D0 + te < 0 && Di <= Dv + te, (Di Inf)/(
  D0 Ip + Inf (Dv + te)),
  D0 >= Dv + te && te < 0 && D0 + te >= 0 && Di <= Dv + te, (Di Inf)/(
  Dv Inf + Inf te - Ip te),
  D0 < Dv + te && D0 + Di >= Dv + te && D0 + te < 0 && 
   Di <= Dv + te, (Inf (D0 + Di - Dv - te))/(D0 (Inf + Ip)),
  D0 < Dv + te && D0 + Di >= Dv + te && te < 0 && D0 + te >= 0 && 
   Di <= Dv + te, (Inf (D0 + Di - Dv - te))/(D0 Inf - Ip te),
  te == 0 && D0 < Dv && D0 + Di >= Dv, (D0 + Di - Dv)/D0,
  D0 < Dv + te && D0 + Di >= Dv + te && te >= 0, 
  1 + (Di - Dv)/(D0 - te),
  D0 >= Dv + te && D0 + te < 0 && Di > Dv + te && Dv + te > 0, (
  Inf (Dv + te))/(D0 Ip + Inf (Dv + te)),
  D0 >= Dv + te && D0 + te >= 0 && Di > Dv + te && Dv + te > 0, (
  Inf (Dv + te))/(Dv Inf + (Inf - Ip) te),
  True, 0]
infectable[t_] := UnitStep[t]*UnitStep[D0 - t]
{$Assumptions = 
   Di < Dv && Dv < D0 && Di > 0 && Di \[Element] Reals && 
    D0 \[Element] Reals && Inf > 0 && Inf \[Element] Reals && 
    Ip > 0  && Ip \[Element] Reals  };
A = FullSimplify[(Inf/N/D0)*phi*tau*
   Integrate[
    Integrate[
     Integrate[
      infectable[ti] Cp[te, Di, Dv, D0], {tx, te + Dv, ti + Di}], {ti,
       te + Dv - Di, te + Dv}], {te, -Dv, D0 - Dv + Di}]]

It used to work when the Cp was defined as (Cor) for the integral (A).

Cor[te_, Di_, Dv_, D0_] :=
 Which[Di > Dv + te && Dv + te > 0, 1,
  te >= 0 && Dv + te <= D0, Di/Dv,
  te < 0 && Di <= Dv + te && Dv + te <= D0, Di/(Dv + te),
  te < 0 && Dv + te > D0, (Di - Dv - te + D0)/D0,
  Di + D0 < Dv + te || Dv + te <= 0, 0,
  True, 1 + (-Di + Dv)/(te - D0)]

infectable[t_] := UnitStep[t]*UnitStep[D0 - t]

{$Assumptions = 
   Di < Dv && Dv < D0 && Di > 0 && Di \[Element] Reals && 
    D0 \[Element] Reals && Inf > 0 && Inf \[Element] Reals && 
    Ip > 0  && Ip \[Element] Reals  };
A = FullSimplify[(Inf/N/D0)*phi*tau*
   Integrate[
    Integrate[
     Integrate[
      infectable[ti] Cor[te, Di, Dv, D0], {tx, te + Dv, 
       ti + Di}], {ti, te + Dv - Di, te + Dv}], {te, -Dv, 
     D0 - Dv + Di}]]

Here is the result from that equation

1/(4 D0 Dv N)
  Inf phi tau (Di (2 D0 Di^2 + (3 Di - 2 Dv) (Di - Dv) Dv) + 
    2 Dv ((-Di + Dv)^3 Log[1 - Di/Dv] + Di^3 Log[Dv/Di]))
Attachments:
POSTED BY: Welling Oei
6 Replies
Sort By:

It used to work when the Cp was defined as (Cor) for the integral (A).

But here the inner integration over tx is trivial because tx is not a variable of the integrand

In[18]:= FreeQ[infectable[ti] Cor[te, Di, Dv, D0], tx]
Out[18]= True

so one could do as well

In[23]:= A1 = 
 FullSimplify[(Inf/N/D0)*phi*tau*
   Integrate[
       Integrate[(ti + Di - te - Dv) infectable[ti] Cor[te, Di, Dv,  D0], {ti, te + Dv - Di, te + Dv}], 
    {te, -Dv, D0 - Dv + Di}]]
Out[23]= (Inf phi tau (Di (2 D0 Di^2 + (3 Di - 2 Dv) (Di - Dv) Dv) + 
   2 Dv ((-Di + Dv)^3 Log[1 - Di/Dv] + Di^3 Log[Dv/Di])))/(4 D0 Dv N)

in the more general case this is also true

In[28]:= FreeQ[infectable[ti] Cp[te, Di, Dv, D0], tx]
Out[28]= True

but the observation does not bring anything. How did you check the consistency (lack of contradictions) of this:

Cp[te_, Di_, Dv_, D0_] := 
 Which[Di + D0 <= Dv + te || Dv + te <= 0,                                      0, 
  D0 < Dv + te && D0 + te < 0 && Di > Dv + te,                                  Inf/(Inf + Ip), 
  D0 < Dv + te && D0 + te >= 0 && Di > Dv + te,                                 (D0 Inf)/(D0 Inf - Ip te), 
  Di <= Dv + te && ((te == 0 && D0 >= Dv) || (D0 >= Dv + te && te >= 0)),       Di/Dv, 
  D0 >= Dv + te && D0 + te < 0 && Di <= Dv + te,                                (Di Inf)/(D0 Ip + Inf (Dv + te)), 
  D0 >= Dv + te && te < 0 && D0 + te >= 0 && Di <= Dv + te,                     (Di Inf)/(Dv Inf + Inf te - Ip te), 
  D0 < Dv + te && D0 + Di >= Dv + te && D0 + te < 0 &&  Di <= Dv + te,          (Inf (D0 + Di - Dv - te))/(D0 (Inf + Ip)), 
  D0 < Dv + te && D0 + Di >= Dv + te && te < 0 && D0 + te >= 0 && Di <= Dv + te,(Inf (D0 + Di - Dv - te))/(D0 Inf - Ip te), 
  te == 0 && D0 < Dv && D0 + Di >= Dv,                                          (D0 + Di - Dv)/D0,  
  D0 < Dv + te && D0 + Di >= Dv + te && te >= 0,                                1 + (Di - Dv)/(D0 - te), 
  D0 >= Dv + te && D0 + te < 0 && Di > Dv + te &&  Dv + te > 0,                 (Inf (Dv + te))/(D0 Ip + Inf (Dv + te)), 
  D0 >= Dv + te && D0 + te >= 0 && Di > Dv + te &&  Dv + te > 0,                (Inf (Dv + te))/(Dv Inf + (Inf - Ip) te), 
  True,                                                                         0
]

if the default evaluates to 0, then the first case is useless, because it also evaluates to 0.

Have you tried to reduce that using your Assumptions:

{$Assumptions = 
   Di < Dv && Dv < D0 && Di > 0 && Di \[Element] Reals && 
    D0 \[Element] Reals && Inf > 0 && Inf \[Element] Reals && Ip > 0 &&
     Ip \[Element] Reals};

there you state Di < Dv < D0 and in Cp you consider te == 0 && D0 < Dv. Di, Dv, and D0 seem to be some constants. te is the dynamic variable. Despite the fact that Inf and Ip remain undefined, plots like

Plot[Cp[x,1,2,4], {x, -1, -5}]

do give a continous function. That suggests that most cases in Cp are somehow useless in the Di < Dv < D0 range and possibly after simplification you could reach a result. If all that Inf and Ip related stuff is kicked out in a mood of jazz it came to

In[56]:= Cp[te_, Di_, Dv_, D0_] := 
 Which[Di <= Dv + te && ((te == 0 && D0 >= Dv) || (D0 >= Dv + te && te >= 0)),  Di/Dv,
  te == 0 && D0 < Dv && D0 + Di >= Dv, (D0 + Di - Dv)/D0, 
  D0 < Dv + te && D0 + Di >= Dv + te && te >= 0, 1 + (Di - Dv)/(D0 - te), 
  True, 0]

In[57]:= infectable[t_] := UnitStep[t]*UnitStep[D0 - t]

In[58]:= {$Assumptions = 
   Di < Dv && Dv < D0 && Di > 0 && Di \[Element] Reals && 
    D0 \[Element] Reals && Inf > 0 && Inf \[Element] Reals && Ip > 0 &&
     Ip \[Element] Reals};

In[59]:= A = 
 FullSimplify[(Inf/N/D0)*phi*tau*
   Integrate[
    Integrate[
     Integrate[
      infectable[ti] Cp[te, Di, Dv, D0], {tx, te + Dv, ti + Di}], {ti,
       te + Dv - Di, te + Dv}], {te, -Dv, D0 - Dv + Di}]]

Out[59]= (Inf phi tau (Di (6 D0 Di^2 + 
      Dv (5 Di^2 - 15 Di Dv + 6 Dv^2)) - 
   12 (Di - Dv)^3 Dv ArcTanh[Di/(Di - 2 Dv)]))/(12 D0 Dv N)
POSTED BY: Udo Krause

I will suggest removing the nested Integrate[Integrate[...]] and simply use one Integrate with all ranges inside. That way Integrate will know about ranges for all variables of integration, not just the one it is handling at the time. This might make a difference.
POSTED BY: Daniel Lichtblau

You're welcome. The specific setting where Do>Dv>Di and te<0

In[7]:= (* First setting: Do>Dv>Di and te<0 *)
Clear[Cp, infectable]
Cp[te_,Di_,Dv_,D0_]:=Which[D0>=Dv+te&&te<0&&D0+te>=0&&Di<=Dv+te,(Di Inf)/(Dv Inf+Inf te-Ip te),D0<Dv+te&&D0+Di>=Dv+te&&te<0&&D0+te>=0&&Di<=Dv+te,(Inf (D0+Di-Dv-te))/(D0 Inf-Ip te),D0>=Dv+te&&D0+te>=0&&Di>Dv+te&&Dv+te>0,(Inf (Dv+te))/(Dv Inf+(Inf-Ip) te),True,0]
infectable[t_]:=UnitStep[t]*UnitStep[D0-t]
In[6]-= {$Assumptions=Di<Dv&&Dv<D0&&Di>0&&Di\[Element]Reals&&D0\[Element]Reals&&Inf>0&&Inf\[Element]Reals&&Ip>0&&Ip\[Element]Reals};
In[10]:= A=FullSimplify[(Inf/N/D0)*phi*tau*Integrate[Integrate[Integrate[infectable[ti] Cp[te,Di,Dv,D0],{tx,te+Dv,ti+Di}],{ti,te+Dv-Di,te+Dv}],{te,-Dv,D0-Dv+Di}]]

has a result (shown as picture because of the formatting )

first setting

because of the factor in the denominator (Inf - Ip)^4 it's clear that a condition appears. That's solved by Mathematica. If one uses one Integrate only, one integral remains.
POSTED BY: Udo Krause

Cp is indeed the complete (not all relevant) conditions for various settings we consider in the calculation. The full formulas for Cp was derived by Mathematica when I specified:

I'm still not confident with the model, for example it could be written equivalently as

In[19]:= Clear[infectable, preinfect, Di, Dv, D0, Inf, Ip]
{$Assumptions = Di < Dv && Di > 0 && D0 > 0};
infectable[t_] := UnitStep[t]*UnitStep[D0 - t]
preinfect[t_] := infectable[-t]
(* NRp=Inf/D0*FullSimplify[Integrate[infectable[ti],{ti,te,te+Dv}]]+\
Ip/D0*FullSimplify[Integrate[preinfect[ti],{ti,te,te+Dv}]] *)
NRp = Integrate[
   Inf infectable[ti] + Ip preinfect[ti], {ti, te, te + Dv}, 
   Assumptions -> te \[Element] Reals]/D0

Out[23]= (1/D0)(Inf UnitStep[
    Dv + te] (UnitStep[
       D0 - te] (D0 - 
        te + (-D0 + Dv + te) UnitStep[D0 - Dv - te]) UnitStep[
       te] - (-1 + 
        UnitStep[
         te]) ((Dv + te) UnitStep[-Dv - te] UnitStep[
          D0 - Dv - 
           te] + (D0 + (-D0 + Dv + te) UnitStep[
             D0 - Dv - te]) UnitStep[Dv + te])) + 
  Ip UnitStep[
    D0 + Dv + 
     te] ((-te + (Dv + te) UnitStep[-Dv - te]) UnitStep[-te] UnitStep[
       D0 + te] - (-1 + 
        UnitStep[
         D0 + te]) ((D0 + Dv + te) UnitStep[-Dv - te] UnitStep[-D0 - 
           Dv - te] + (D0 + (Dv + te) UnitStep[-Dv - te]) UnitStep[
          D0 + Dv + te])))

if this NRp is inverted, one gets a clear statement where the NRp is zero (there 1/NRp equals to ComplexInfinity) which was absent in other formulations ...

1/NRp reformulated
POSTED BY: Udo Krause

Dear Udo, Thank you for your kind assessment and suggestions.

Cp is indeed the complete (not all relevant) conditions for various settings we consider in the calculation. The full formulas for Cp was derived by Mathematica when I specified:

{$Assumptions = Di < Dv && Di > 0 && D0 > 0};
infectable[t_] := UnitStep[t]*UnitStep[D0 - t]
preinfect[t_] := UnitStep[D0 + t]*UnitStep[-t]
NRp = Inf/D0*
    FullSimplify[Integrate[infectable[ti], {ti, te, te + Dv}]] + 
   Ip/D0*FullSimplify[Integrate[preinfect[ti], {ti, te, te + Dv}]];
FullSimplify[PiecewiseExpand[TR/NRp], 
 Assumptions -> 
  Di < Dv && Di > 0 && D0 > 0 && Inf > 0 && Ip > 0  &&  D0 > 0  && 
   te \[Element] Reals]

Following your suggestion, I tried to reduce the conditions to a specific setting where Do>Dv>Di and te<0:

Cp[te_, Di_, Dv_, D0_] := 
 Which[D0 >= Dv + te && te < 0 && D0 + te >= 0 && 
   Di <= Dv + te, (Di Inf)/(Dv Inf + Inf te - Ip te), 
  D0 < Dv + te && D0 + Di >= Dv + te && te < 0 && D0 + te >= 0 && 
   Di <= Dv + te, (Inf (D0 + Di - Dv - te))/(D0 Inf - Ip te), 
  D0 >= Dv + te && D0 + te >= 0 && Di > Dv + te && Dv + te > 0, (
  Inf (Dv + te))/(Dv Inf + (Inf - Ip) te), True, 0]
infectable[t_] := UnitStep[t]*UnitStep[D0 - t]

{$Assumptions = 
   Di < Dv && Dv < D0 && Di > 0 && Di \[Element] Reals && 
    D0 \[Element] Reals && Inf > 0 && Inf \[Element] Reals && Ip > 0 &&
     Ip \[Element] Reals};

A = FullSimplify[(Inf/N/D0)*phi*tau*
   Integrate[
    Integrate[
     Integrate[
      infectable[ti] Cp[te, Di, Dv, D0], {tx, te + Dv, ti + Di}], {ti,
       te + Dv - Di, te + Dv}], {te, -Dv, D0 - Dv + Di}]]

Or for the setting where Do>Di>Dv and te<0:

Cp[te_, Di_, Dv_, D0_] := 
 Which[D0 >= Dv + te && te < 0 && D0 + te >= 0 && 
   Di <= Dv + te, (Di Inf)/(Dv Inf + Inf te - Ip te), 
  D0 < Dv + te && D0 + Di >= Dv + te && te < 0 && D0 + te >= 0 && 
   Di <= Dv + te, (Inf (D0 + Di - Dv - te))/(D0 Inf - Ip te), 
  D0 >= Dv + te && D0 + te >= 0 && Di > Dv + te && Dv + te > 0, (
  Inf (Dv + te))/(Dv Inf + (Inf - Ip) te), True, 0]
infectable[t_] := UnitStep[t]*UnitStep[D0 - t]

{$Assumptions = 
   Di > Dv && Dv < D0 && Di > 0 && Di \[Element] Reals && 
    D0 \[Element] Reals && Inf > 0 && Inf \[Element] Reals && Ip > 0 &&
     Ip \[Element] Reals};

A = FullSimplify[(Inf/N/D0)*phi*tau*
   Integrate[
    Integrate[
     Integrate[
      infectable[ti] Cp[te, Di, Dv, D0], {tx, te + Dv, ti + Di}], {ti,
       te + Dv - Di, te + Dv}], {te, -Dv, D0 - Dv + Di}]]

But I still got a conditional expression as the output.

Do you have some suggestions how to solve this?

Kind regards, Welling
POSTED BY: Welling Oei

In the second case for the setting where Do > Di > Dv and te < 0 the conditional expression under the remaining integral works into your hands stating

-Inf (Dv + te)^2 (-2 Di + Dv + te)/(2 (Dv Inf + (Inf-Ip) te)) for D0>= Dv + te && Di > Dv + te && Dv + te > 0
0 elsewhere

because you say te < 0 the upper integration limitD0 +Di -Dv for te must be smaller than 0

second setting

there is a contradiction: If Di > Dv and D0 > Dv then D0 + Di - Dv > 0 so te > 0 (not te < 0) on the upper limit .... okay, let's send this back to you and continue with the integration limits:

The lower limit is fullfilled: D0 > 0 && Di > 0 && 0 >= 0 because D0> Dv in the assumptions and if te < 0 at least on the lower limit one has Dv > 0.

The upper limit D0 >= D0 + Di && Di > D0 + Di && D0 + Di > 0 is not fullfilled because Di > 0; one has to determine where the integrand is 0 and where not. To get out of this mess one could plot it

With[{Dv = 1, Di = 4.5, D0 = 8},
 Plot[If[D0 >= Dv + te && Di > Dv + te && Dv + te > 0, 
   1, -1], {te, -Dv, D0 + Di + Dv}]
 ]

fortunately there is only one switch from 1 to -1 at te = Di - Dv leaving us behind with the integral

In[30]:= Clear[Dv, Di, D0]
Inf phi tau Integrate[-((Inf (Dv + te)^2 (-2 Di + Dv + te))/(
    2 (Dv Inf + (Inf - Ip) te))), {te, -Dv, Di - Dv}, 
   Assumptions -> Dv > 0]/(D0 N)

Out[31]= ConditionalExpression[-((
  Inf^2 phi tau (-Di (Inf - Ip) (4 Di^2 (Inf - Ip)^2 - 6 Dv^2 Ip^2 + 
        9 Di Dv Ip (-Inf + Ip)) + 
     6 Dv^2 Ip^2 (2 Di (Inf - Ip) + Dv Ip) Log[(Dv Ip)/(
       Di Inf - Di Ip + Dv Ip)]))/(12 D0 (Inf - Ip)^4 N)), 
 Di Ip^2 < Di Inf Ip || Di (Inf - Ip) (Di (Inf - Ip) + Dv Ip) < 0]

the integrals are done. Please check for errors.
POSTED BY: Udo Krause
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract