Mercurial > repos > public > sbplib

changeset 512:4ef2d2a493f1 feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
add H_xyz to hypsyst to compute H-norm
author Ylva Rydin <ylva.rydin@telia.com>
date Mon, 26 Jun 2017 19:23:19 +0200
parents 57f3493f851b
children bc39bb984d88
files +scheme/Hypsyst3dCurve.m +time/+expint/Magnus_4.m +time/+expint/Magnus_mp.m
diffstat 3 files changed, 9 insertions(+), 3 deletions(-) [+]
line wrap: on
line diff
--- a/+scheme/Hypsyst3dCurve.m	Thu Jun 08 10:33:36 2017 +0200
+++ b/+scheme/Hypsyst3dCurve.m	Mon Jun 26 19:23:19 2017 +0200
@@ -216,6 +216,8 @@
             obj.Hxii = kr(I_n, ops_xi.HI, I_eta,I_zeta);
             obj.Hetai = kr(I_n, I_xi, ops_eta.HI,I_zeta);
             obj.Hzetai = kr(I_n, I_xi,I_eta, ops_zeta.HI);
+            Hi = obj.Hxii*obj.Hetai*obj.Hzetai;
+            obj.H = inv(Hi);
             
             obj.index_w = (kr(ops_xi.e_l, O_eta,O_zeta)==1);
             obj.index_e = (kr(ops_xi.e_r, O_eta,O_zeta)==1);
--- a/+time/+expint/Magnus_4.m	Thu Jun 08 10:33:36 2017 +0200
+++ b/+time/+expint/Magnus_4.m	Mon Jun 26 19:23:19 2017 +0200
@@ -11,7 +11,9 @@
    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;
+  % v = expm(Omega)*v;
+     toler = 10^(-8);
+  v = time.expint.expm_Arnoldi(-Omega,v,k,toler,100);
 else
    %v = krylov(k*D,v);
    v = expm(k*D)*v;
--- a/+time/+expint/Magnus_mp.m	Thu Jun 08 10:33:36 2017 +0200
+++ b/+time/+expint/Magnus_mp.m	Mon Jun 26 19:23:19 2017 +0200
@@ -5,10 +5,12 @@
 
 if isa(D,'function_handle')
    % v = krylov(k*D(t +k/2*t),v);
-   v = expm(k*D(t +k/2))*v;
+%   v = expm(k*D(t +k/2))*v;
+      toler = 10^(-5);
+   expm_Arnoldi(-D,v,k,toler,100)
 else
    %v = krylov(k*D,v);
-   v = expm(k*D)*v;
+   % v = expm(k*D)*v;
 end
 
 end
\ No newline at end of file