3
$\begingroup$

I am trying to get the integral of a rather complicated function: $$ \int_{t_1}^{t_2}m_{\text{basal}}+m_{\text{size}}\exp\left(\frac{1}{2}\left(\ln(S_{t_1})+\frac{u-t_1}{t_2-t_1}(\ln(S_{t_2})-\ln(S_{t_1}))+\frac{\delta^2(t_2-u)(u-t_1)}{t_2-t_1}\right)\right) \mathrm{d} u $$ But since I had the feeling that it might be solvable by a computer algebra system, I plugged it into Mathematica and got a solution. However, the solution involves a square root of negative numbers and is therefore complex, which surprised me. The integrand is a real function and evaluation of the solution as well, since all complex parts cancel out at the end, magically. But I would like to avoid any complex parts. As a mathematica newbie, I am not too familiar with its functionalities. Is there any way to enforce the symbolic integration just using real functions? Or putting in some constraints of variables, like $t1<t2$, in order not to get a square-root of a negative number?

Even more problematic is that evaluating the solution showed complete nonsense results, when compared to numerical integration of the Integrand. The results of the numeric integration is in a range, I would have expected (due to values of the integrand). Any ideas why this happens?

Below you can find the code:

betaExample := 0.1
deltaExample := 0.1
mbasalExample := 0.1
msizeExample := 0.1
S0Example := 1.0
t1Example := 0.0
t2Example := 1.0
S1Example := 1.105170918
S2Example := 4.0552

exampleSub = {beta -> betaExample, delta -> deltaExample, 
  msize -> msizeExample, mbasal -> mbasalExample, t1 -> t1Example, 
  t2 -> t2Example, S1 -> S1Example, S2 -> S2Example};

examplifyExpression[expr_] := expr /. exampleSub
Integrand = mbasal + msize* Exp[(1/2)*(Log[S1] + ((u - t1)/(t2 - t1))*(Log[S2] - Log[S1]) + 
                                (delta^2/4)*((t2 - u)*(u - t1))/(t2 - t1))]
(* mbasal + E^(1/2 ((delta^2 (t2 - u) (-t1 + u))/(4 (-t1 + t2)) + Log[S1] + 
            ((-t1 + u) (-Log[S1] + Log[S2]))/(-t1 + t2))) msize *)

AnaInt1 = FullSimplify[Integrate[Integrand, {u, t1, t2}]]
(* mbasal (-t1 + t2) + (2 Sqrt[2] msize Sqrt[t1 - t2] (-Sqrt[S2] 
DawsonF[(delta^2 (t1 - t2) - 4 Log[S1] + 4 Log[S2])/(4 Sqrt[2] delta Sqrt[t1 - t2])] + 
Sqrt[S1]
DawsonF[(delta^2 (-t1 + t2) - 4 Log[S1] + 4 Log[S2])/(4 Sqrt[2] delta Sqrt[t1 - t2])]))/delta *)

examplifyExpression[AnaInt1]
(* 6.57826*10^23 + 0. I *)

NIntegrate[examplifyExpression[Integrand], {u, t1Example, t2Example}]
(* 0.206923 *)
$\endgroup$
1
  • $\begingroup$ Welcome to the Mathematica Stack Exchange. Please present a minimal example that can be copied/pasted/executed in our own notebook sessions. $\endgroup$
    – Syed
    Jan 22 at 13:23

2 Answers 2

4
$\begingroup$

You have precision issues. Use arbitrary precision, for example:

AnaInt1 /. exampleSub
(* 0.1 + 0. I *)

AnaInt1 /. SetPrecision[exampleSub, 10]
(* 0.2481047 *)

The result now matches the NIntegrate one:

NIntegrate[Integrand /. exampleSub, {u, t1Example, t2Example}]
(* 0.248105 *)

Regarding the output not being real: The integral Mathematica returned is real for $t_2\geq t_1$, despite the presence of terms like $\sqrt{t_1-t_2}$. Think, for example, of a simpler function: $$ f(t_1,t_2) = \sqrt{t_1-t_2}\left(\sqrt{t_1-t_2-1} + \sqrt{t_1-t_2-2}\right).$$ Again, it may look this function is not real for $t_2\geq t_1$, but it is! $$ f(t_1,t_2) = i\sqrt{t_2-t_1}\left(i\sqrt{t_2-t_1+1} + i\sqrt{t_2-t_1+2}\right) = -\sqrt{t_2-t_1}\left(\sqrt{t_2-t_1+1} + \sqrt{t_2-t_1+2}\right).$$

The same goes for your function: $i$ coming from sqare root will cancel with the one(s) coming from Erfi.

What you can do to get a different form, is to use Refine:

Refine[AnaInt1, {S1, S2, delta} ∈ PositiveReals && t2 >= t1]

You can cancel the $i$s by noting that:

FunctionExpand[DawsonF[x I]]
(* 1/2 I E^x^2 Sqrt[π] Erf[x] *)
$\endgroup$
2
  • $\begingroup$ Thanks for your answer. Was not aware of these kind of issues. I am not only interested in the evaluation of this but also in the analytical expression itself. Do you also know an answer to the first part? Like one can force a solution with real numbers, by giving constraints on the inputs or so? $\endgroup$ Jan 22 at 15:09
  • $\begingroup$ @VincentWieland, please see my update. $\endgroup$
    – Domen
    Jan 22 at 15:33
2
$\begingroup$
$Version

(* "14.0.0 for Mac OS X ARM (64-bit) (December 13, 2023)" *)

Clear["Global`*"]

For complicated expressions, avoid machine precision numbers.

betaExample = 1/10;
deltaExample = 1/10;
mbasalExample = 1/10;
msizeExample = 1/10;
S0Example = 1;
t1Example = 0;
t2Example = 1;
S1Example = 1.105170918 // Rationalize[#, 0] &;
S2Example = 4.0552 // Rationalize[#, 0] &;

exampleSub = {beta -> betaExample, delta -> deltaExample, 
   msize -> msizeExample, mbasal -> mbasalExample, t1 -> t1Example, 
   t2 -> t2Example, S1 -> S1Example, S2 -> S2Example};

examplifyExpression[expr_] := expr /. exampleSub

Integrand = 
  mbasal + msize*
     Exp[(1/2)*(Log[
          S1] + ((u - t1)/(t2 - t1))*(Log[S2] - Log[S1]) + (delta^2/
            4)*((t2 - u)*(u - t1))/(t2 - t1))] // Simplify;

Using FullSimplify with the following is too slow for me.

AnaInt1 = Simplify[Integrate[Integrand, {u, t1, t2}]]

(* -mbasal t1 + mbasal t2 + (1/delta)
 E^(-((delta^4 (t1 - t2)^2 + 16 Log[S1]^2 - 32 Log[S1] Log[S2] + 
    16 Log[S2]^2)/(32 delta^2 (t1 - t2)))) msize Sqrt[2 π] S1^(1/4) S2^(
  1/4) Sqrt[
  t1 - t2] (-Erfi[(delta^2 (t1 - t2) - 4 Log[S1] + 4 Log[S2])/(
      4 Sqrt[2] delta Sqrt[t1 - t2])] + 
    Erfi[(delta^2 (-t1 + t2) - 4 Log[S1] + 4 Log[S2])/(
     4 Sqrt[2] delta Sqrt[t1 - t2])]) *)

res1 = examplifyExpression[AnaInt1] // FullSimplify

(* 1/10 + 1/5 (789975451153/70506939)^(1/4) E^(
  1/3200 + 50 Log[357399673791/97402773125]^2)
   Sqrt[π] (Erf[(1 - 400 Log[357399673791/97402773125])/(40 Sqrt[2])] + 
    Erf[(1 + 400 Log[357399673791/97402773125])/(40 Sqrt[2])]) *)

res1 // N[#, 15] &

(* 0.248104731870805 *)

res2 = NIntegrate[examplifyExpression[Integrand], {u, t1Example, t2Example}, 
  WorkingPrecision -> 15]

(* 0.248104731870805 *)

res1 - res2

(* 0.*10^-16 *)
$\endgroup$
1
  • $\begingroup$ Yeah, it also took hours on my laptop. $\endgroup$ Jan 22 at 15:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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