Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger1dCurve.m @ 495:b91d23271481 feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Added new penalty parameter to the interface in shrodinger curve
author Ylva Rydin <ylva.rydin@telia.com>
date Fri, 24 Feb 2017 09:04:02 +0100
parents 6b8297f66c91
children 437fad4a47e1
comparison
equal deleted inserted replaced
494:c6d03f951a7f 495:b91d23271481
22 gamm 22 gamm
23 end 23 end
24 24
25 methods 25 methods
26 % Solving SE in the form u_t = i*u_xx +i*V on deforming 1D domain; 26 % Solving SE in the form u_t = i*u_xx +i*V on deforming 1D domain;
27 function obj = Schrodinger1dCurve(m,order,V,constJi) 27 function obj = Schrodinger1dCurve(g,m,order,V,constJi)
28 default_arg('V',0); 28 default_arg('V',0);
29 default_arg('constJi',false) 29 default_arg('constJi',false)
30 xilim={0 1}; 30 xilim={0 1};
31 if constJi 31 if constJi
32 ops = sbp.D2Standard(m,xilim,order); 32 ops = sbp.D2Standard(m,xilim,order);
43 obj.M = ops.M; 43 obj.M = ops.M;
44 obj.e_l = ops.e_l; 44 obj.e_l = ops.e_l;
45 obj.e_r = ops.e_r; 45 obj.e_r = ops.e_r;
46 obj.d1_l = ops.d1_l; 46 obj.d1_l = ops.d1_l;
47 obj.d1_r = ops.d1_r; 47 obj.d1_r = ops.d1_r;
48 obj.grid = g;
48 49
49 50
50 if isa(V,'function_handle') 51 if isa(V,'function_handle')
51 V_vec = V(obj.x); 52 V_vec = V(obj.x);
52 else 53 else
105 end 106 end
106 107
107 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 108 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
108 % u denotes the solution in the own domain 109 % u denotes the solution in the own domain
109 % v denotes the solution in the neighbour domain 110 % v denotes the solution in the neighbour domain
110 [e_u,d_u,s_u] = obj.get_boundary_ops(boundary); 111 [e_u,d_u,s_u,p_u] = obj.get_boundary_ops(boundary);
111 [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 112 [e_v,d_v,s_v,p_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
112 113
113 a = -s_u* 1/2 * 1i ; 114 a1 = -s_u* 1/2 * 1i ;
114 b = a'; 115 b1 = a1';
116 gamma = -1*s*a(p_u,p_u)/2;
115 117
116 tau = b*d_u; 118 tau = b1*d_u;
117 sig = -a*e_u; 119 sig = -a1*e_u;
118 120
119 closure = obj.Hi * (tau*e_u' + sig*d_u'); 121 closure = @(a) obj.Hi * (tau*e_u' + sig*d_u' + gamma(a));
120 penalty = obj.Hi * (-tau*e_v' - sig*d_v'); 122 penalty = @(a) obj.Hi * (-tau*e_v' - sig*d_v' - gamma(a));
121 end 123 end
122 124
123 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 125 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
124 % The right boundary is considered the positive boundary 126 % The right boundary is considered the positive boundary
125 function [e,d,s,p] = get_boundary_ops(obj,boundary) 127 function [e,d,s,p] = get_boundary_ops(obj,boundary)