Mercurial > repos > public > sbplib

comparison +time/+expint/Magnus_4.m @ 513:bc39bb984d88 feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added arnoldi krylov subspace approximation
author Ylva Rydin <ylva.rydin@telia.com>
date Mon, 26 Jun 2017 20:15:54 +0200
parents 4ef2d2a493f1
children
comparison
equal deleted inserted replaced
512:4ef2d2a493f1 513:bc39bb984d88
1 % Takes one time step of size k using a fourth order magnus integrator 1 % Takes one time step of size k using a fourth order magnus integrator
2 % starting from v_0 and where the function F(v,t) gives the 2 % starting from v_0 and where the function F(v,t) gives the
3 % time derivatives. 3 % time derivatives.
4 function v = Magnus_4(v,D, t , k) 4 function v = Magnus_4(v, D, t , k , matrixexp ,tol)
5
6
5 7
6 if isa(D,'function_handle') 8 if isa(D,'function_handle')
7 % v = krylov(k*D(t +k/2*t),v); 9 c1 = 1/2 - sqrt(3)/6;
8 c1 = 1/2 - sqrt(3)/6; 10 c2 = 1/2 + sqrt(3)/6;
9 c2 = 1/2 + sqrt(3)/6; 11
10 12 A1 = D(t +c1*k);
11 A1 = D(t +c1*k); 13 A2 = D(t + c2*k);
12 A2 = D(t + c2*k); 14 Omega = 1/2*(A1 + A2) + sqrt(3)*k/12*(A1*A2-A2*A1);
13 Omega = k/2*(A1 + A2) + sqrt(3)*k^2/12*(A1*A2-A2*A1);
14 % v = expm(Omega)*v;
15 toler = 10^(-8);
16 v = time.expint.expm_Arnoldi(-Omega,v,k,toler,100);
17 else 15 else
18 %v = krylov(k*D,v); 16 Omega = D;
19 v = expm(k*D)*v;
20 end 17 end
21 18
19
20 switch matrixexp
21 case 'expm'
22 v = expm(k*Omega)*v;
23 case 'Arnoldi'
24 v = time.expint.expm_Arnoldi(-Omega,v,k,tol,100);
25 otherwise
26 error('No such matrix exponential evaluation')
27
22 end 28 end
29 end