Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger1dCurve.m @ 506:2047ff9640b3 feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
tried to make a symmetric scheme
author Ylva Rydin <ylva.rydin@telia.com>
date Mon, 29 May 2017 13:54:16 +0200
parents b7280c6c6b9a
children 1341e8580d64
comparison
equal deleted inserted replaced
505:49bc2929dfb1 506:2047ff9640b3
15 ddt_x_r 15 ddt_x_r
16 ddt_x_l 16 ddt_x_l
17 a 17 a
18 a_xi 18 a_xi
19 Ji 19 Ji
20 J
20 t_up 21 t_up
21 x 22 x
22 23
23 V_mat 24 V_mat
24 D1 25 D1
83 % neighbour_scheme is an instance of Scheme that should be interfaced to. 84 % neighbour_scheme is an instance of Scheme that should be interfaced to.
84 % neighbour_boundary is a string specifying which boundary to interface to. 85 % neighbour_boundary is a string specifying which boundary to interface to.
85 86
86 function [D] = d_fun(obj,t) 87 function [D] = d_fun(obj,t)
87 obj.variable_update(t); % In driscretization? 88 obj.variable_update(t); % In driscretization?
88 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); 89 D = sqrt(obj.Ji)*(-0.5*(obj.D1*obj.a*sqrt(obj.Ji) + obj.a*obj.D1*sqrt(obj.Ji)) + 1i*obj.D2(diag(obj.Ji))*obj.Ji + 1i*obj.V_mat);
89 90 % 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);
90 end 91 end
91 92
92 93
93 function [] = variable_update(obj,t) 94 function [] = variable_update(obj,t)
94 if (t == obj.t_up) 95 if (t == obj.t_up)
100 ddt_x_l = obj.ddt_x_l(t); 101 ddt_x_l = obj.ddt_x_l(t);
101 obj.x = obj.xi*(x_r -x_l) + x_l; 102 obj.x = obj.xi*(x_r -x_l) + x_l;
102 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))); 103 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)));
103 104
104 obj.Ji = sparse(diag(1./(x_r - x_l + 0*obj.x))); 105 obj.Ji = sparse(diag(1./(x_r - x_l + 0*obj.x)));
106 obj.J = sparse(x_r -x_l);
105 obj.a_xi = sparse(diag(-1*(ddt_x_r - ddt_x_l + 0*obj.x))); 107 obj.a_xi = sparse(diag(-1*(ddt_x_r - ddt_x_l + 0*obj.x)));
106 obj.t_up = t; 108 obj.t_up = t;
107 end 109 end
108 end 110 end
109 111
114 [e,d,s,p] = obj.get_boundary_ops(boundary); 116 [e,d,s,p] = obj.get_boundary_ops(boundary);
115 117
116 switch type 118 switch type
117 % Dirichlet boundary condition 119 % Dirichlet boundary condition
118 case {'D','d','dirichlet'} 120 case {'D','d','dirichlet'}
119 tau1 = @(t) s * 1i*obj.Ji(p,p)^2*d; 121 tau1 = @(t) s * 1i*obj.Ji(p,p)*d;
120 tau2 = @(t) obj.Ji*(-1*s*obj.a(p,p) - abs(obj.a(p,p)))/4*e; 122 tau2 = @(t) (1*s*obj.a(p,p))/2*e;
121 closure = @(t) obj.Hi*tau1(t)*e' + obj.Hi*tau2(obj.a)*e'; 123 closure = @(t)(obj.Ji*obj.Hi*tau1(t)*e' + obj.Hi*tau2(obj.a)*e');
122 124
123 switch class(data) 125 switch class(data)
124 case 'double' 126 case 'double'
125 penalty = @(t) -obj.Ji*(obj.Hi*tau1*data+obj.Hi*tau2(obj.a)*data); 127 penalty = @(t) -(obj.Hi*tau1*data+obj.Hi*tau2(obj.a)*data);
126 % case 'function_handle' 128 % case 'function_handle'
127 % penalty = @(t)-obj.Hi*tau*data(t); 129 % penalty = @(t)-obj.Hi*tau*data(t);
128 otherwise 130 otherwise
129 error('Wierd data argument!') 131 error('Weird data argument!')
130 end 132 end
131 133
132 % Unknown, boundary condition 134 % Unknown, boundary condition
133 otherwise 135 otherwise
134 error('No such boundary condition: type = %s',type); 136 error('No such boundary condition: type = %s',type);