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 *)
Improve this question
$\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$ Commented Jan 22, 2024 at 13:23

2 Answers 2

Reset to default
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] *)
Improve this answer
$\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$ Commented Jan 22, 2024 at 15:09
  • $\begingroup$ @VincentWieland, please see my update. $\endgroup$ Commented Jan 22, 2024 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 *)
Improve this answer
$\endgroup$
1
  • $\begingroup$ Yeah, it also took hours on my laptop. $\endgroup$ Commented Jan 22, 2024 at 15:09

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.