# NIntegration of trace of matrix, which was obtained by NDSolve, converging too slowly

e.g.

ClearAll["Global*"]

H[a1_, a2_, a3_, t_] := {{0, Cos[a1*t], Cos[a2*t]}, {Cos[a1*t], 0, Cos[a3*t]}, {Cos[a2*t], Cos[a3*t], 0}}

Ut[a1_, a2_, a3_] := Ut[a1, a2, a3] = NDSolveValue[{Derivative[1][u][x] == (-I)*H[a1, a2, a3, x] . u[x], u[0] == IdentityMatrix[3]}, u, {x, 0, 2*Pi}]

S[t_] := Piecewise[{{CosIntegral[20.*2.*Pi*Abs[t]] - CosIntegral[0.0001*2.*Pi*Abs[t]], t != 0}, {0, t == 0}}]

Mi[a1_, a2_, a3_] := (1/2)*NIntegrate[Tr[Ut[a1, a2, a3][2*Pi] . ConjugateTranspose[Ut[a1, a2, a3][t1]] . {{0, 1/3, 1/3}, {1/3, 0, 1/3}, {1/3, 1/3, 0}} . Ut[a1, a2, a3][t1] .
ConjugateTranspose[Ut[a1, a2, a3][t2]] . {{0, 1/3, 1/3}, {1/3, 0, 1/3}, {1/3, 1/3, 0}} . Ut[a1, a2, a3][t2] . {{0, 0, 0}, {0, 1/Sqrt[2], 0}, {0, 0, -(1/Sqrt[2])}}]*S[t1 - t2],
{t1, 0, 2*Pi}, {t2, 0, t1}]


but Mi[1,2,3] warning Numerical integration converging too slowly.

Could it be modify?

• Function ut[1,2,3]doesn't evaluate. Perhaps the initial conditions should be u[0] == {1, 1, 1} ? Dec 21 '21 at 12:40
• @Ulrich Neumann MatrixForm@H should be removed. Sorry for my input error, question have been edited. Dec 21 '21 at 12:50
• Perhaps Method -> "InterpolationPointsSubdivision" inside NIntegrate helps. Dec 21 '21 at 12:58
• @Ulrich Neumann it outcome 3.78277 - 1.2844 I for Mi[1,2,3], but warning NIntegrate obtained 7.56554 -2.5688 I and 0.000363099 for the integral and error estimates. Dec 22 '21 at 1:36
• In addition to what has been suggested already, it's worth trying Ut = ParametricNDSolveValue[{Derivative[1][u][x] == (-I)* H[a1, a2, a3, x] . u[x], u[0] == IdentityMatrix[3]}, u, {x, 0, 2*Pi}, {a1, a2, a3}]. ParametricNDSolveValue can do derivatives w.r.t. the parameters and should speed up the numerical solutions of the differential equation. Dec 23 '21 at 8:51

The definition of S[t] is numerically problematic because of the condition t==0 and t != 0.

Try a smoothed version S[t]  which is continuous near t==0:

S = Function[{t}, Which[t == 0, Log[20/10^-4], True,CosIntegral[20*2*Pi*t] - CosIntegral[10^-4*2*Pi*t] ] // Evaluate]


Evaluation time of the integral can be decreased significantly (factor3) if you integrate inside a region:

reg=DiscretizeRegion[ImplicitRegion[{0 <= t1 <= 2 Pi && 0 <= t2 <= t1}, {t1, t2}]]

Mi[a1_, a2_, a3_] := (1/2)*
NIntegrate[
Tr[Ut[a1, a2, a3][2*Pi] .
ConjugateTranspose[
Ut[a1, a2, a3][t1]] . {{0, 1/3, 1/3}, {1/3, 0, 1/3}, {1/3, 1/3,
0}} . Ut[a1, a2, a3][t1] .
ConjugateTranspose[
Ut[a1, a2, a3][t2]] . {{0, 1/3, 1/3}, {1/3, 0, 1/3}, {1/3, 1/3,
0}} . Ut[a1, a2, a3][
t2] . {{0, 0, 0}, {0, 1/Sqrt[2], 0}, {0, 0, -(1/Sqrt[2])}}]*
S[t1 - t2], Element[{t1, t2}, reg],
Method -> "InterpolationPointsSubdivision"]

Mi[1,2,3]
(*3.78283 - 1.28444 I*)


Further improvement (evaluation speed) might be achieved if the integrand is interpolated inside reg

ipMii[a1_, a2_, a3_, n_ (*number of interpolation points*)] :=
Block[{t1, t2, ip,
reg = ImplicitRegion[{0 <= t1 <= t2 && 0 <= t2 <= t1}, {t1, t2}]},
ip = FunctionInterpolation[(1/
2) (Tr[Ut[a1, a2, a3][2*Pi] .
ConjugateTranspose[
Ut[a1, a2, a3][t1]] . {{0, 1/3, 1/3}, {1/3, 0, 1/3}, {1/3,
1/3, 0}} . Ut[a1, a2, a3][t1] .
ConjugateTranspose[
Ut[a1, a2, a3][t2]] . {{0, 1/3, 1/3}, {1/3, 0, 1/3}, {1/3,
1/3, 0}} .
Ut[a1, a2, a3][
t2] . {{0, 0, 0}, {0, 1/Sqrt[2], 0}, {0,
0, -(1/Sqrt[2])}}]*S[t1 - t2]), {t1, 0, 2 Pi}, {t2, 0,
2 Pi}, InterpolationPoints -> n, InterpolationOrder -> 1]] ;
NIntegrate[ipMii[1, 2, 3, 20][t1, t2], Element[{t1, t2}, reg],Method -> "InterpolationPointsSubdivision"]
(*3.86797 - 1.33425 I*)
`
• WOW! Thank you @Ulrich Neumann for your detailed answer! That's really a lot can learned. Dec 23 '21 at 1:36
• Mi wroks in my laptop. When I run Mii, however, MMA warns The integrand has evaluated to non-numerical values for all sampling points. I have tried NIntegrate[ip[t1, t2], {t1,0,Pi}, {t2,0,Pi}]] and it works. It would be nice if you can tell me what was wrong when I copied the codes.@Ulrich Neumann Dec 23 '21 at 7:48
• I didn't understand the messages too, but the result seems to be ok. Dec 23 '21 at 8:09
• See my modified answer. Dec 23 '21 at 8:20
• Thumbs-up！Thank you again, @Ulrich Neumann. It is an excellent solution. Dec 23 '21 at 8:48