# Implementing the Peano-Baker Series with Mathematica

I would like to calculate the Peano-Baker series with mathematica, namely an integral expansion to approximate the state transition matrix:

Say I have a matrix $$A(t)$$ and I want to calculate:

$$\Phi(t,t_0) = e^{\int_{t_0}^t A(\sigma)d \sigma}$$

To approximately calculate this I can expand the Peano-Baker Series as follows:

$$\Phi(t,t_0) = I + \int_{t_0}^t A(\sigma_1)d\sigma_1+ \int_{t_0}^{t} A(\sigma_1) \int_{t_0}^t A(\sigma_2) d \sigma_2 \sigma_1+ \ldots \int_{t_0}^{t} A(\sigma_1) \int_{t_0}^t A(\sigma_2) \dots \int_{t_0}^{\sigma_{j-1}}A(\sigma_j)d\sigma_j \dots d\sigma_2 d \sigma_1 + \ldots$$

What would be a smart way to do this with Mathematica? I would like a function of this sort:

PeanoBakerSeries[A_?MatrixQ,t_,t0_,j_?IntegerQ]:=Module[{n, dummies, PB},
n = Length@A;
dummies = Table[Symbol["\[Sigma]" <> ToString[i]], {i, 1, k}];
....]


I was also wondering if there is an undocumented function for this (I bet there is!).

• MatrixExp[Integrate[A[s],{s, t0, t}]]? – Henrik Schumacher Aug 1 '19 at 12:05
• Your form for the PB series looks different from the one I am accustomed to (e.g. the one here). You can use @Henrik's expression if the appropriate commutator is zero, but in general, you will want to exploit Cauchy's iterated integral formula. – J. M. will be back soon Aug 1 '19 at 12:15

I actually feel like I figured out a possible solution:

     PeanoBakerSeries[A_,t_,t0_,k_?IntegerQ]:=Module[{n,dummies,PB},n=Length@A[t];
dummies=Table[Symbol["\[Sigma]"<>ToString[i]],{i,1,k}];
PB=Integrate[A,{dummies[],t0,t}];
PB=IdentityMatrix[n];
Do[With[{intervals=FlattenAt[{{dummies[],t0,t},Sequence@@@Table[{{dummies[[i]],t0,dummies[[i-1]]}},{i,2,j}]},2]},
PB=PB+Integrate[Product[testA[dummies[[i]]],{i,1,j}],Sequence@@intervals]],{j,1,k}];
PB
]


Nevertheless I will keep the question open, in case someone feels to suggest something smarter or an undocumented function that does this.