0
$\begingroup$

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?

$\endgroup$
6
  • $\begingroup$ Function ut[1,2,3]doesn't evaluate. Perhaps the initial conditions should be u[0] == {1, 1, 1} ? $\endgroup$ Dec 21 '21 at 12:40
  • $\begingroup$ @Ulrich Neumann MatrixForm@H should be removed. Sorry for my input error, question have been edited. $\endgroup$ Dec 21 '21 at 12:50
  • $\begingroup$ Perhaps Method -> "InterpolationPointsSubdivision" inside NIntegrate helps. $\endgroup$ Dec 21 '21 at 12:58
  • $\begingroup$ @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. $\endgroup$ Dec 22 '21 at 1:36
  • $\begingroup$ 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. $\endgroup$ Dec 23 '21 at 8:51
0
$\begingroup$

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*)
$\endgroup$
6
  • $\begingroup$ WOW! Thank you @Ulrich Neumann for your detailed answer! That's really a lot can learned. $\endgroup$ Dec 23 '21 at 1:36
  • $\begingroup$ 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 $\endgroup$ Dec 23 '21 at 7:48
  • $\begingroup$ I didn't understand the messages too, but the result seems to be ok. $\endgroup$ Dec 23 '21 at 8:09
  • $\begingroup$ See my modified answer. $\endgroup$ Dec 23 '21 at 8:20
  • $\begingroup$ Thumbs-up!Thank you again, @Ulrich Neumann. It is an excellent solution. $\endgroup$ Dec 23 '21 at 8:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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