Mercurial > repos > public > sbplib

changeset 508:508b7493be94 feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
moved all multiplications and divisions bu J and Ji to the siscretization to easier find the error
author Ylva Rydin <ylva.rydin@telia.com>
date Fri, 02 Jun 2017 11:33:04 +0200
parents 1341e8580d64
children 566304a5733e
files +scheme/Schrodinger1dCurve.m
diffstat 1 files changed, 11 insertions(+), 8 deletions(-) [+]
line wrap: on
line diff
--- a/+scheme/Schrodinger1dCurve.m	Mon May 29 14:07:55 2017 +0200
+++ b/+scheme/Schrodinger1dCurve.m	Fri Jun 02 11:33:04 2017 +0200
@@ -17,6 +17,7 @@
         a
         a_xi
         Ji
+        J
         t_up
         x
         
@@ -85,8 +86,9 @@
         
         function [D] = d_fun(obj,t)
             obj.variable_update(t); % In driscretization?
-            D= obj.Ji*(-0.5*(obj.D1*obj.a - obj.a_xi + obj.a*obj.D1) + 1i*obj.D2(diag(obj.Ji)) + 1i*obj.V_mat);
-            
+             D = (-0.5*(obj.D1*obj.a + obj.a*obj.D1) + 1i*obj.D2(diag(obj.Ji)) + 1i*obj.V_mat);
+%             D = (-0.5*(obj.D1*obj.a -obj.a_xi+ obj.a*obj.D1) + 1i*obj.D2(diag(obj.Ji)) + 1i*obj.V_mat);
+          %   D= obj.Ji*(-sqrt(obj.a)*obj.D1*sqrt(obj.a) + 0.5*obj.a_xi + 1i*obj.D2(diag(obj.Ji)) + 1i*obj.V_mat);
         end
         
         
@@ -102,6 +104,7 @@
                 obj.a = sparse(diag((-ddt_x_l*( x_r - x_l) - (obj.x-x_l)*(ddt_x_r-ddt_x_l))/(x_r-x_l)));            
                 
                 obj.Ji = sparse(diag(1./(x_r - x_l + 0*obj.x)));
+                obj.J = sparse(x_r -x_l);
                 obj.a_xi = sparse(diag(-1*(ddt_x_r - ddt_x_l + 0*obj.x)));
                 obj.t_up = t;
             end
@@ -116,17 +119,17 @@
             switch type
                 % Dirichlet boundary condition
                 case {'D','d','dirichlet'}
-                    tau1 = @(t) s * 1i*obj.Ji(p,p)^2*d;
-                    tau2 = @(t) obj.Ji*(-1*s*obj.a(p,p) - abs(obj.a(p,p)))/4*e;
-                    closure = @(t) obj.Hi*tau1(t)*e' + obj.Hi*tau2(obj.a)*e';
+                    tau1 = @(t) s * 1i*obj.Ji(p,p)*d;
+                    tau2 = @(t) (1*s*obj.a(p,p))/2*e;
+                    closure = @(t)(obj.Hi*tau1(t)*e' + obj.Hi*tau2(obj.a)*e');
                     
                     switch class(data)
                         case 'double'
-                            penalty = @(t) -obj.Ji*(obj.Hi*tau1*data+obj.Hi*tau2(obj.a)*data);
+                            penalty = @(t) -(obj.Hi*tau1*data+obj.Hi*tau2(obj.a)*data);
                             %                      case 'function_handle'
                             %                           penalty = @(t)-obj.Hi*tau*data(t);
                         otherwise
-                            error('Wierd data argument!')
+                            error('Weird data argument!')
                     end
                     
                     % Unknown, boundary condition
@@ -143,7 +146,7 @@
             
             a1 =   s_u* 1/2 * 1i ;
             b1 =   -s_u* 1/2 * 1i;
-            gamma = @(a) -obj.Ji*s_u*a(p_u,p_u)/2*e_u;
+            gamma = @(a) -s_u*a(p_u,p_u)/2*e_u;
             
             tau = @(t) a1*obj.Ji(p_u,p_u)^2*d_u;
             sig = b1*e_u;