How to make integrals of interpolated matrices constructed by other interpolated quantities?

Question

I have the following code producing an interpolated unitary matrix:

H[t_] := PauliMatrix[3] + Cos[t]*PauliMatrix[1];
sol = First[NDSolve[{I*D[matU[s],s] == H[s].matU[s],matU[0] == IdentityMatrix[2]}, matU, {s, 0, 20}]];
mU = matU /. sol;

This produces the unitary matrix mU which I can use as a function. For example, writing mU[4] gives a matrix. However, if I want to integrate this, for example writing NIntegrate[mU[t], {t, 1, 5}] it does not work. It says

NIntegrate: Integrand InterpolatingFunction[{{0.,20.}},{5,3,1,{282},{4},0,0,0,0,Automatic,<<3>>},{{<<1>>}},{{{{1. +0. I,0. +0. I},{0. +0. I,1. +0. I}},{{0. -1. I,0. -1. I},{0. -1. I,0. +1. I}}},<<9>>,<<272>>},{Automatic}][t] is not numerical at {t} = {1.00015}.

The first question is: how do I do this integral?

Second part: I use the unitary matrix to define time-dependent vectors as follows,

vec = {0, 1};
vect[t_] := mU[t].vec;
ρ[t_] := Outer[Times, vect[t], vect[t]\[Conjugate]];

Then, I try to integrate using NIntegrate[ρ[t], {t, 1, 5}] and I get two different errors:

Outer: Heads Conjugate and Dot at positions 3 and 2 are expected to be the same.

NIntegrate: Integrand Outer[Times,<<1>>,Conjugate[InterpolatingFunction[{{0.,20.}},{5,3,1,{282},{4},0,0,0,0,Automatic,<<3>>},{{0.,<<9>>,<<272>>}},{{{{<<2>>},{<<2>>}},{{<<2>>},{<<2>>}}},{{{<<2>>},{<<2>>}},{{<<2>>},{<<2>>}}},<<8>>,<<272>>},{Automatic}][t].{0,1}]] is not numerical at {t} = {1.00015}.

Why do I get these errors and how can I avoid them?

Thanks everybody for the help!

Edit1: When Integrate is directly applied to the original interpolated matrix it seems to work. Therefore, the value of NIntegrate[ρ[t], {t, 1, 5}] can be obtained as follows. First, we make Mathematica numerically solve the differential equation for ρ:

solρ = First[NDSolve[{I*D[ρS[s], s] == 
      H[s].ρS[s] - ρS[s].H[s], ρS[0] == 
      Outer[Times, vec, vec]}, ρS, {s, 0, 20}]];
ρt = ρS /.solρ;

Then, we can integrate it using Integrate[ρt[t], {t, 1, 5}]. However, if we try to compute Integrate[Sin[t]*ρt[t], {t, 1, 5}] or even Integrate[Sin[t]+ρt[t], {t, 1, 5}] we do not get any number as output (but we also do not get error messages).

The same numerical result can be obtained with the following code which follows the one given in an answer:

vect2[t_, i_] := mU[t][[i, 2]];
ρ2[t_, i_, j_] := vect2[t, i]*vect2[t, j]\[Conjugate];
Table[NIntegrate[ρ2[t, i, j], {t, 1, 5}], {i, 2}, {j, 2}]

When using this, Mathematica does complain a bit but arrives at the same result up to 10^(-7). Moreover, with this code we can calculate the following Table[NIntegrate[Sin[t]*ρ2[t, i, j], {t, 1, 5}], {i, 2}, {j, 2}]. However, this code seems quite a workaround, which I would prefer to avoid in more complex situations.

Answer 1

There seems to be a problem with NIntegrate and matrices (although for simple matrices NIntegrate works.).

You may try calculating the different matrix elements separately:

H[t_] := PauliMatrix[3] + Cos[t]*PauliMatrix[1];
sol = First[
   NDSolve[{I*D[matU[s], s] == H[s] . matU[s], 
     matU[0] == IdentityMatrix[2]}, matU, {s, 0, 20}]];
mU = matU /. sol;
vect[t_] := mU[t][[ All,2]];
Clear[[Rho]];
[Rho][t_?NumericQ, i_, j_] := vect[t][[i]] Conjugate[vect[t]][[j]];

With this we can calculate the matrix:

Table[NIntegrate[\[Rho][t, i, j], {t, 1, 5}] , {i, 2}, {j, 2}]

Answer 2

NIntegrate is expecting the integrand to evaluate to a number, not a matrix, which is why it doesn't work.

Instead of using NIntegrate, why don't you just augment your original differential equation? For example:

{int, m} = NDSolveValue[
    {
        I matU'[s] == H[s] . matU[s], matU[0] == IdentityMatrix[2],
        i'[s] == matU[s], i[0] == {{0, 0}, {0, 0}}
    },
    {i, matU},
    {s, 0, 20}
];

Then:

int[5] - int[1]

{{-0.00620054 + 0.650422 I, -0.925243 - 0.481496 I}, {0.925243 - 0.481496 I, -0.00620054 - 0.650422 I}}

which agrees with:

Integrate[m[s], {s, 1, 5}]

{{-0.00620039 + 0.650422 I, -0.925243 - 0.481496 I}, {0.925243 - 0.481496 I, -0.00620039 - 0.650422 I}}

Similarly, for your ρ example, you can do:

{int, sinint, ρ} = NDSolveValue[
    {
        I ρS'[s] == H[s] . ρS[s] - ρS[s] . H[s], ρS[0] == Outer[Times, vec, vec],
        i'[s] == ρS[s], i[0] == {{0, 0}, {0, 0}},
        isin'[s] == Sin[s] ρS[s], isin[0] == {{0, 0}, {0, 0}}
    },
    {i, isin, ρS},
    {s, 0, 20}
];

Then:

int[5] - int[1]

{{1.43151 + 0. I, 0.147931 + 0.453175 I}, {0.147931 - 0.453175 I, 2.56849 + 0. I}}

which agrees with Daniel's answer, while:

sinint[5] - sinint[1]

{{0.499208 + 0. I, -0.231473 + 0.138672 I}, {-0.231473 - 0.138672 I, -0.242568 + 0. I}}

gives the value of:

Integrate[Sin[t] ρ[t], {t, 1, 5}]

Note that using NDSolveValue is much more convenient than using NDSolve.