# Disagreement between symbolic and numerical integration

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, 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 *)

• Welcome to the Mathematica Stack Exchange. Please present a minimal example that can be copied/pasted/executed in our own notebook sessions.
– Syed
Jan 22 at 13:23

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] *)

• 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? Jan 22 at 15:09
• @VincentWieland, please see my update. Jan 22 at 15:33
\$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 *)
`
• Yeah, it also took hours on my laptop. Jan 22 at 15:09