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