# NIntegrate on a solution of a matrix ODE

I've seen similar questions on this site but somehow the solutions there didn't manage to solve my specific problem.

I have a function mat1 that takes a square $$n \times n$$ matrix G, and some final time tfinal, and solves the following ODE numerically: $$u'(t) = G(t) u(t)$$ $$u(0) = \mathrm{id}_{n\times n}$$ The code is:

mat1[G_, tfinal_] := Block[{t}, NDSolveValue[{u'[t] == G[t].u[t], u[0] == IdentityMatrix[Dimensions[G[0]][[1]]]}, u, {t, 0, tfinal},
Method -> "ExplicitRungeKutta"]]


Let's take an example matrix-valued function $$g(t)$$:

g[t_?NumericQ] := {{Sin[t], 0}, {Cos[t], t}}


Mathematica has no problems solving the ODE with g as the input matrix:

mat1[g, 10][1.21]
(*Result: {{1.90977, 0.}, {1.92296, 2.07912}}*)


But when I want to numerically integrate it, I get the following error:

NIntegrate[mat1[g, 10][t], {t, 0, 10}]
(*NIntegrate::inum: Integrand InterpolatingFunction[{{0.,10.}},{5,3,1,{98},{4},0,0,0,0,Automatic,{},{},False},{{0.,0.120666,0.60333,0.874901,<<43>>,6.97746,7.05172,7.12517,<<48>>}},{{{{1.,0.},{0.,1.}},{{0.,0.},{1.,0.}}},{{{1.0073,0.},{0.121253,1.00731}},{{0.121252,0.},{1.0146,0.121548}}},<<48>>,<<48>>},{Automatic}][t] is not numerical at {t} = {0.000960178}.*)
(*NIntegrate::inum: Integrand InterpolatingFunction[{{0.,10.}},{5,3,1,{98},{4},0,0,0,0,Automatic,{},{},False},{{0.,0.120666,0.60333,0.874901,<<43>>,6.97746,7.05172,7.12517,<<48>>}},{{{{1.,0.},{0.,1.}},{{0.,0.},{1.,0.}}},{{{1.0073,0.},{0.121253,1.00731}},{{0.121252,0.},{1.0146,0.121548}}},<<48>>,<<48>>},{Automatic}][t] is not numerical at {t} = {0.000960178}.*)


I've also tried defining a function in between:

mat2[t_?NumericQ] := mat1[g, 10][t]


But I get the same error:

NIntegrate[mat2[t], {t, 0, 10}]
(*NIntegrate::inum: Integrand mat2[t] is not numerical at {t} = {0.0795732}.*)


It looks like even with the NumericQ, Mathematica is trying to manipulate the integrand with a symbolic $$t$$ before putting numbers in.

EDIT:

It looks like the above code works fine for a real-valued function, as opposed to matrices:

mat1[G_, tfinal_] := Block[{t}, NDSolveValue[{u'[t] == G[t]*u[t], u[0] ==1},u, {t, 0, tfinal}, Method -> "ExplicitRungeKutta"]]

g[t_?NumericQ] := Sin[t]

NIntegrate[mat1[g, 10][t], {t, 0, 10}]

(*Result: 36.4662*)


So it looks like the problem has something to do with $$g$$ being a matrix. I'm not sure how though.

• Welcome to Mathematica.SE, Sahand! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – Chris K Apr 21 '19 at 19:48
• This is an answer to a related question. – Anton Antonov Apr 22 '19 at 7:04

As the error message says, the problem is that mat2[0.0795732] is not numerical. It is instead a 2x2 matrix of numbers. You could do something like:

mat2[t_?NumericQ] := mat1[g, 10][t][[1,1]]
NIntegrate[mat2[t], {t, 0, 10}]


36.4662

On the other hand, it is much simpler to just have NDSolveValue do the integration for you:

mat1[G_,tfinal_] := NDSolveValue[
{
int'[t] == u[t], int[0] == ConstantArray[0, Dimensions[G[0]]],
u'[t]==G[t].u[t], u[0]==IdentityMatrix[Dimensions[G[0]][[1]]]
},
{u, int},
{t,0,tfinal}
]


Then:

mat1[g, 10]


and:

mat1[g, 10][[2]][10]


{{36.4662, 0.}, {3.69638*10^20, 5.23821*10^20}}

agreeing with the above result.

• I see! So Mathematica also counts matrices with numerical elements as "non-numerical". So I have to integrate it element-wise, or use NDSolve. So is there no direct way to integrate matrices with NIntegrate itself? – Sahand Tabatabaei Apr 21 '19 at 18:40
• @SahandTabatabaei Yes, there is a direct way to integrate matrices with NIntegrate. I will post a related answer after a day or two. (I plan to extend the rule described here or program new one...) – Anton Antonov Apr 21 '19 at 21:51

You can use Integrate to directly antidifferentiate an interpolating function. If $$f(t)$$ is an interpolating function with domain $$(a,b)$$, Integrate[f[t], t] returns an interpolating function with the same domain equal to $$\int_a^t f(\tau) \; d\tau\,.$$

To get the definite integral, plug the end point:

Integrate[mat1[g, 10][t], t] /. t -> 10
(*  {{36.4662, 0.}, {3.69611*10^20, 5.23781*10^20}}  *)

• You could also use Derivative[-1][mat1[g, 10]] to construct the interpolating function without an argument, e.g., Derivative[-1][mat1[g, 10]][10] so that no ReplaceAll is needed. – Carl Woll Apr 21 '19 at 19:50
• @CarlWoll Thanks. I was going to add that if the OP wanted greater accuracy, using InterpolationOrder -> Allin NDSolve would likely produce a more accurate integral (by either of our methods, I suppose), but that option does not work with matrix ODEs and certain methods, such "ExplicitRungeKutta" with a difference order greater than 3. (Just reported as [CASE:4249898].) – Michael E2 Apr 21 '19 at 20:09