Mercurial > repos > public > sbplib
changeset 495:b91d23271481 feature/quantumTriangles
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 | c6d03f951a7f |
| children | 437fad4a47e1 |
| files | +scheme/Schrodinger.m +scheme/Schrodinger1dCurve.m |
| diffstat | 2 files changed, 19 insertions(+), 14 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/Schrodinger.m Thu Feb 23 09:32:30 2017 +0100 +++ b/+scheme/Schrodinger.m Fri Feb 24 09:04:02 2017 +0100 @@ -22,12 +22,15 @@ methods % Solving SE in the form u_t = i*u_xx -i*V; - function obj = Schrodinger(m,xlim,order,V) + function obj = Schrodinger(g,order,V) default_arg('V',0); - ops = sbp.D2Standard(m,xlim,order); + obj.grid = g; + m = N(obj.grid); + obj.x = points(obj.grid); + ops = sbp.D2Standard(m,{obj.x(1) obj.x(end)},order); - obj.x=ops.x; - obj.h=ops.h; + + obj.h = ops.h; obj.D2 = ops.D2; obj.H = ops.H; obj.Hi = ops.HI;
--- a/+scheme/Schrodinger1dCurve.m Thu Feb 23 09:32:30 2017 +0100 +++ b/+scheme/Schrodinger1dCurve.m Fri Feb 24 09:04:02 2017 +0100 @@ -24,7 +24,7 @@ methods % Solving SE in the form u_t = i*u_xx +i*V on deforming 1D domain; - function obj = Schrodinger1dCurve(m,order,V,constJi) + function obj = Schrodinger1dCurve(g,m,order,V,constJi) default_arg('V',0); default_arg('constJi',false) xilim={0 1}; @@ -45,6 +45,7 @@ obj.e_r = ops.e_r; obj.d1_l = ops.d1_l; obj.d1_r = ops.d1_r; + obj.grid = g; if isa(V,'function_handle') @@ -107,17 +108,18 @@ function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) % u denotes the solution in the own domain % v denotes the solution in the neighbour domain - [e_u,d_u,s_u] = obj.get_boundary_ops(boundary); - [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); - - a = -s_u* 1/2 * 1i ; - b = a'; + [e_u,d_u,s_u,p_u] = obj.get_boundary_ops(boundary); + [e_v,d_v,s_v,p_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); - tau = b*d_u; - sig = -a*e_u; + a1 = -s_u* 1/2 * 1i ; + b1 = a1'; + gamma = -1*s*a(p_u,p_u)/2; - closure = obj.Hi * (tau*e_u' + sig*d_u'); - penalty = obj.Hi * (-tau*e_v' - sig*d_v'); + tau = b1*d_u; + sig = -a1*e_u; + + closure = @(a) obj.Hi * (tau*e_u' + sig*d_u' + gamma(a)); + penalty = @(a) obj.Hi * (-tau*e_v' - sig*d_v' - gamma(a)); end % Ruturns the boundary ops and sign for the boundary specified by the string boundary.