Mercurial > repos > public > sbplib
changeset 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 | 32a24485f3e8 |
| files | +time/+expint/Magnus_4.m +time/+expint/Magnus_mp.m +time/+expint/Magnus_mp.m~ +time/+expint/expm_Arnoldi.m +time/Magnus4.m +time/MagnusMP.m |
| diffstat | 6 files changed, 119 insertions(+), 24 deletions(-) [+] |
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
--- a/+time/+expint/Magnus_mp.m Mon Jun 26 19:23:19 2017 +0200 +++ b/+time/+expint/Magnus_mp.m Mon Jun 26 20:15:54 2017 +0200 @@ -1,16 +1,21 @@ % Takes one time step of size k using the magnus midpoinr % starting from v_0 and where the function F(v,t) gives the % time derivatives. -function v = Magnus_mp(v,D, t , k) +function v = Magnus_mp(v,D, t , k,matrixexp,tol) if isa(D,'function_handle') - % v = krylov(k*D(t +k/2*t),v); -% v = expm(k*D(t +k/2))*v; - toler = 10^(-5); - expm_Arnoldi(-D,v,k,toler,100) + Omega = D(t +k/2); 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+expint/Magnus_mp.m~ Mon Jun 26 20:15:54 2017 +0200 @@ -0,0 +1,17 @@ +% Takes one time step of size k using the magnus midpoinr +% starting from v_0 and where the function F(v,t) gives the +% time derivatives. +function v = Magnus_mp(v,D, t , k,matrixexp,tol) + +if isa(D,'function_handle') + switch matrixexp + case 'expm' + v = expm(k*D(t +k/2))*v; + case 'Arnol' + v = time.expint.expm_Arnoldi(-D(t +k/2),v,k,toler,100); +else + %v = krylov(k*D,v); + % v = expm(k*D)*v; +end + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+expint/expm_Arnoldi.m Mon Jun 26 20:15:54 2017 +0200 @@ -0,0 +1,54 @@ +function y = expm_Arnoldi(A,v,t,toler,m) +% +% y = expm_Arnoldi(A,v,t,toler,m) +% +% computes $y = \exp(-t A) v$ +% input: A (n x n)-matrix, v n-vector, t>0 time interval, +% toler>0 tolerance, m maximal Krylov dimension +% +% Copyright (c) 2012 by M.A. Botchev +% Permission to copy all or part of this work is granted, +% provided that the copies are not made or distributed +% for resale, and that the copyright notice and this +% notice are retained. +% +% THIS WORK IS PROVIDED ON AN "AS IS" BASIS. THE AUTHOR +% PROVIDES NO WARRANTY WHATSOEVER, EITHER EXPRESSED OR IMPLIED, +% REGARDING THE WORK, INCLUDING WARRANTIES WITH RESPECT TO ITS +% MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE. +% +n = size (v,1); +V = zeros(n ,m+1); +H = zeros(m+1,m); + +beta = norm(v); +V(:,1) = v/beta; +resnorm = inf; + +j=0; + +while resnorm > toler + j = j+1; + w = A*V(:,j); + for i=1:j + H(i,j) = w'*V(:,i); + w = w - H(i,j)*V(:,i); + end + H(j+1,j) = norm(w); + e1 = zeros(j,1); e1(1) = 1; + ej = zeros(j,1); ej(j) = 1; + s = [0.01, 1/3, 2/3, 1]*t; + for q=1:length(s) + u = expm(-s(q)*H(1:j,1:j))*e1; + beta_j(q) = -H(j+1,j)* (ej'*u); + end + resnorm = norm(beta_j,'inf'); + % fprintf('j = %d, resnorm = %.2e\n',j,resnorm); + if resnorm<=toler + break + elseif j==m + disp('warning: no convergence within m steps'); + end + V(:,j+1) = w/H(j+1,j); +end +y = V(:,1:j)*(beta*u);
--- a/+time/Magnus4.m Mon Jun 26 19:23:19 2017 +0200 +++ b/+time/Magnus4.m Mon Jun 26 20:15:54 2017 +0200 @@ -8,17 +8,23 @@ v m n + matrixexp + tol end methods - function obj = Magnus4(D, k, t0, v0) + function obj = Magnus4(D, k, t0, v0,matrixexp,tol) + default_arg('matrixexp','expm') + default_arg('tol',1e-6) obj.D = D; obj.k = k; obj.t = t0; obj.v = v0; obj.m = length(v0); obj.n = 0; + obj.matrixexp = matrixexp; + obj.tol = tol; end function [v,t] = getV(obj) @@ -27,7 +33,7 @@ end function obj = step(obj) - obj.v = time.expint.Magnus_4(obj.v,obj.D, obj.t, obj.k); + obj.v = time.expint.Magnus_4(obj.v,obj.D, obj.t, obj.k, obj.matrixexp, obj.tol); obj.t = obj.t + obj.k; obj.n = obj.n + 1; end
--- a/+time/MagnusMP.m Mon Jun 26 19:23:19 2017 +0200 +++ b/+time/MagnusMP.m Mon Jun 26 20:15:54 2017 +0200 @@ -8,17 +8,23 @@ v m n + matrixexp + tol end methods - function obj = MagnusMP(D, k, t0, v0) + function obj = MagnusMP(D, k ,t0,v0, matrixexp,tol) + default_arg('matrixexp','expm') + default_arg('tol',1e-6) obj.D = D; obj.k = k; obj.t = t0; obj.v = v0; obj.m = length(v0); obj.n = 0; + obj.matrixexp = matrixexp; + obj.tol = tol; end function [v,t] = getV(obj) @@ -27,7 +33,7 @@ end function obj = step(obj) - obj.v = time.expint.Magnus_mp(obj.v,obj.D, obj.t, obj.k); + obj.v = time.expint.Magnus_mp(obj.v,obj.D, obj.t, obj.k,obj.matrixexp,obj.tol); obj.t = obj.t + obj.k; obj.n = obj.n + 1; end