Mercurial > repos > public > sbplib
diff +time/+expint/Magnus_4.m @ 513:bc39bb984d88 feature/quantumTriangles
Added arnoldi krylov subspace approximation
author | Ylva Rydin <ylva.rydin@telia.com> |
---|---|
date | Mon, 26 Jun 2017 20:15:54 +0200 |
parents | 4ef2d2a493f1 |
children |
line wrap: on
line diff
--- a/+time/+expint/Magnus_4.m Mon Jun 26 19:23:19 2017 +0200 +++ b/+time/+expint/Magnus_4.m Mon Jun 26 20:15:54 2017 +0200 @@ -1,22 +1,29 @@ % Takes one time step of size k using a fourth order magnus integrator % starting from v_0 and where the function F(v,t) gives the % time derivatives. -function v = Magnus_4(v,D, t , k) +function v = Magnus_4(v, D, t , k , matrixexp ,tol) + + if isa(D,'function_handle') - % v = krylov(k*D(t +k/2*t),v); - c1 = 1/2 - sqrt(3)/6; - c2 = 1/2 + sqrt(3)/6; - - A1 = D(t +c1*k); - A2 = D(t + c2*k); - Omega = k/2*(A1 + A2) + sqrt(3)*k^2/12*(A1*A2-A2*A1); - % v = expm(Omega)*v; - toler = 10^(-8); - v = time.expint.expm_Arnoldi(-Omega,v,k,toler,100); + c1 = 1/2 - sqrt(3)/6; + c2 = 1/2 + sqrt(3)/6; + + A1 = D(t +c1*k); + A2 = D(t + c2*k); + Omega = 1/2*(A1 + A2) + sqrt(3)*k/12*(A1*A2-A2*A1); else - %v = krylov(k*D,v); - v = expm(k*D)*v; + Omega = D; end + +switch matrixexp + case 'expm' + v = expm(k*Omega)*v; + case 'Arnoldi' + v = time.expint.expm_Arnoldi(-Omega,v,k,tol,100); + otherwise + error('No such matrix exponential evaluation') + +end end \ No newline at end of file