Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger1dCurve.m @ 507:1341e8580d64 feature/quantumTriangles

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