Calculation the integral of the exponential of a matrix

10 visualizzazioni (ultimi 30 giorni)
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
  4 Commenti
Torsten
Torsten il 20 Lug 2021
E.g. if Mmat is nonsingular, the integral is (expm(t*Mmat) - I)*inv(Mmat)
raha ahmadi
raha ahmadi il 21 Lug 2021
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

Accedi per commentare.

Risposta accettata

David Goodmanson
David Goodmanson il 20 Lug 2021
Modificato: David Goodmanson il 20 Lug 2021
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
  3 Commenti
David Goodmanson
David Goodmanson il 21 Lug 2021
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
raha ahmadi il 22 Lug 2021
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:)

Accedi per commentare.

Più risposte (0)

Categorie

Scopri di più su Linear Algebra in Help Center e File Exchange

Prodotti


Release

R2018b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by