Calculation the integral of the exponential of a matrix

I want to integrate of an exponential of a matrix ( E= exp (Mmat.*(t0-t))). First of all ‘expm’ built-in function in MATLAB has complexity of order N^3 ; while my matrix is a circulant matrix, so I can use a theorem instead that makes complexity of the order of NlogN;

The algorithm is :

( in this algorithm, E=exp(A) which A is a matrix )

By this algorithm I can easily calculate the exp of a matrix but my problem now is the calculation of the integral of this function I mean

load filem.mat
t0=1e-4;
Mmatt=@(t)Mmat.*(t0-t);
col1=Mmatt(:,1);
FFtcol1=fft(col1);
expFF=exp(FFtcol1);
expMt=ifft(expFF);% the first coulmn of exp (mt)
expMmat(:,1)=expMt;
for x=2:322
expMmat(:,x) = circshift(expMDelt,x-1,1);
end
fun=expMmat*Dmat;
f=integral(fun,0,1e-4,'ArrayValued',true);

I don’t know how can I define it as a function. I can use handle function or function defining in a separate file code but how can I do it. Also I can use integral built in finction in matlab.

Also I think I have an another option and calculate the amount of the function in each time and use the trapz function for integration but I think Integral is more suitable for this situation. I attach my code and the data. I really appreciate any help

thanks in advance

[Torsten]

Did you google

integral of a matrix exponential

?

[raha ahmadi]

Thank you for your attention. Yes I did. Also I read the refrence paper which was used by MATLAB. My matrix is a circulant matrix which makes problem easier. In fact I did not want use unitary transform for calculatin exp (diag a))

[Torsten]

E.g. if Mmat is nonsingular, the integral is (expm(t*Mmat) - I)*inv(Mmat)

[raha ahmadi]

Dear Toresten thank you for your responce. My matrix is usually sparse and inv function is really inefficient. Altough I use the other algorithms but I get inappropriate answers With best regards

Answer 1

[David Goodmanson]

Hi raha,

I think you can go with the following technique, which does integration, differentiation etc. at the eigenvalue level. I was casual about it and created the circulant matrices by concatentation which of course would have to be improved for large matrices.

t = 4;
m1 = 2*rand(5,1)-1
M = circ(m1)
expM = expm(M*t)

% calculate N = expm(M*t)
c1 = M(:,1);
lam = fft(c1);
lamfun = exp(lam*t);
i1 = ifft(lamfun);
N = circ(i1)
max(abs(N-expM),[],'all')   % check should be small

% calculate D = d/dt expm(M*t)
c1 = M(:,1);
lam = fft(c1);
lamfun = (exp(lam*t)).*lam;
i1 = ifft(lamfun);
D = circ(i1)
max(abs(D - M*expM),[],'all')
max(abs(D - expM*M),[],'all')

% calculate I = Int expm(M*t) dt   (indefinite integral)
c1 = M(:,1);
lam = fft(c1);
lamfun = (exp(lam*t))./lam;
i1 = ifft(lamfun);
I = circ(i1)
max(abs(I*M -expM),[],'all')
max(abs(M*I -expM),[],'all')

function M = circ(m1)
% create circulant matrix by circular shift of columns
n = size(m1,1);
M = m1;
for k = 1:n-1
   M = [M circshift(m1,k,1)];
end
end

[raha ahmadi]

Dear Goodmanson Your method is very innovative, fast and effisient. I really enjoyed and learned very much. Thank you very much With best wishes

[David Goodmanson]

You're welcome, raha. I would like to thank you too for introducing me, and hopefully others as well, to this way of exponentiating circulant matrices. i was not aware of it, and I was able to do one of my favorite things, turning it into an mfile to save in a folder.

[raha ahmadi]

I 'm very happy to hear that. The reference for the algorithm is this <https://backend.orbit.dtu.dk/ws/portalfiles/portal/5265681/thesis.pdf > I agree with you the algorithm is very helpful. I hope you the best and of course thanks again:)