# Numerical integration of numerical function

I am trying to solve numerically this integral $$F(s)=\int_{t1}^{s}e^{-B|s-t|}\rho_{12}(t)dt$$ where $$t1\leq s\leq t2$$ and $$\rho_{12}(t)$$ is the element of density matrix $$\rho(t)$$ which can be obtained numerically by solving the von Neuman equation $$\dot{\rho}(t)=-i[H(t),\rho(t)]$$ for

$$H_{0}(t)=\begin{pmatrix}h1+v t & -{\it i}g \\{\it i}g & -h1-v t \\\end{pmatrix}.$$

and $$t1$$ and $$t2$$ is initial and final times.

It is easy to obtain the density matrix $$\rho_{12}$$ element by using the code below

        g = 1;
h1 = -50;
h2 = 50;
v = 1;
B = 1;
t1 = 0;
t2 = (h2-h1)/v;

H0[t_] = {{h1+v*t , -I*g},{I*g, -h1-v*t }};

Sol=NDSolveValue[{D[rho[t], t] == -I*(H0[t].rho[t] - rho[t].H0[t]),
rho[t1] == {{1, 0}, {0, 0}}}, rho, {t, t1, t2}];

(*rho12[t]=Sol[t][[1,2]]*)

F[s_]=Integrate[Exp[B*Abs[s-t]]*Sol[t][[1,2]],{t,t1,s}]


But my code doesn't work to calculate the F(s). I was wondering if you would be able to help me. Thank you

• You did not write the code for $H_0(t)$ and did not write the definition for $\rho(t)$ using Mathematica code. You just showed the Latex version for these. It will be better to include the code for these also so one does not to write them and make mistake. Sep 19 at 10:49
• @Nasser Thank you for your comment. During the revision the matrix $H_0(t)$ was removed. Now my code is complete and works. Sep 19 at 11:06
• @RadmehrJafari There is a typo in your code with B usage. In Latex equation you use $e^{-B|t2-t|}$, while in the code Exp[B*Abs[t2-t]]. With Exp[-B*Abs[t2-t]] the answer is 0.00034735 + 0.00734009 I. Sep 19 at 12:15
• @Alex Trounev Thank you so much for your comments. I have corrected my question and latex Sep 19 at 12:39
• Changing the question after getting answers is discouraged. It tends to make the current answers seem irrelevant. It's better to ask new question, copying relevant code. A small addition might be ok, but the original question should remain. See mathematica.meta.stackexchange.com/questions/1808/… for more. Sep 19 at 12:48

One "book" method (that is, as recommended in the documentation) is to use Indexed instead Part and the "InterpolationPointsSubdivision" method of NIntegrate:

PrintTemporary@Dynamic@{Clock[Infinity]};
NIntegrate[Exp[B*Abs[t2 - t]]*Indexed[Sol[t], {1, 2}], {t, t1, t2},
Method -> {"InterpolationPointsSubdivision",
"SymbolicProcessing" -> 0,
"MaxSubregions" -> Length@Sol@"Grid"}] // AbsoluteTiming

(*  {252.614, 2.744*10^39 + 2.74278*10^41 I}  *)


Slice out the {1, 2} component and reinterpolate: much faster, if only because we're computing one component instead of four at each function evaluation.

sol12 = Interpolation[Transpose@{
Sol@"Grid",
Sol["ValuesOnGrid"][[All, 1, 2]],
Sol'["ValuesOnGrid"][[All, 1, 2]]}];

PrintTemporary@Dynamic@{Clock[Infinity]};
NIntegrate[Exp[B*Abs[t2 - t]]*sol12[t], {t, t1, t2},
Method -> {"InterpolationPointsSubdivision",
"MaxSubregions" -> Length@ sol12@"Grid"}] // AbsoluteTiming

(*  {49.1007, 2.744*10^39 + 2.74278*10^41 I}  *)


Note that the length of the interpolation "Grid" is nearly 24000, at each node of which there is a weak singularity that slows convergence.

Also, while solving matrix ODEs is convenient for setup, there is a cost deferred in having a more expensive function to evaluate, especially if one wants to deal with the components separately. In a case like this, the advantage of @Nasser's approach should be evident.

Here's a way to separate all the component into interpolating functions (sol is a 2x2 array):

sol = Map[
Prepend,
{#,
Sol["Grid"]}] &,
Transpose[{
Sol["ValuesOnGrid"],
Sol'["ValuesOnGrid"]},
{4, 3, 1, 2}],
{2}];

• Thank you so much for your comments. My question now is that how I can calculate the function F(s). Sep 19 at 12:37

If I understand you right, you could always start easy. No need to make everything very short. Later, once it is working, you can always improve and shorten the code if needed. So you could do

g = 1;
h1 = -50;
h2 = 50;
v = 1;
B = 1;
t1 = 0;
t2 = (h2 - h1)/v
H0[t_] := {{h1 + v*t, -I*g}, {I*g, -h1 - v*t}}
rho = {{rho11[t], rho12[t]}, {rho21[t], rho22[t]}}
(ode = Thread[Flatten /@ (D[rho, t] == -I*(H0[t] . rho - rho . H0[t]))]) // Column


These are the 4 odes to solve. Now setup the IC and call the solver

ic = {rho11[t1] == 1, rho12[t1] == 0, rho21[t1] == 0, rho22[t1] == 0}
Sol = NDSolveValue[{ode, ic}, {rho11, rho12, rho21, rho22}, {t, t1, t2}]


Now integrate $$\rho_{1,2}$$ which is the second one above

 NIntegrate[Exp[B*Abs[t2 - t]]*Sol[[2]][t], {t, t1, t2}]


ps. if the sign meant to be negative then the answer is

NIntegrate[Exp[-B*Abs[t2 - t]]*Sol[[2]][t], {t, t1, t2}]


I do not know if this makes sense or not, you know the Physics. Mathematica also gives warning about convergence.

• Thank you so much Sep 19 at 11:30
• Module[{t}, NIntegrate[Exp[B*Abs[t2 - t]]*Sol[t][[1, 2]], {t, t1, t2}]] gives the same result 3.50298*10^39 + 2.75472*10^41 I directly! Sep 19 at 11:47
• @Ulrich Neumann Thank you for your comment. I have corrected my question. Sep 19 at 13:57