Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger2dCurve.m @ 700:de7f665e1d00 feature/optim

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
add hamiltonian closures
author Ylva Rydin <ylva.rydin@telia.com>
date Mon, 16 Oct 2017 16:51:22 +0200
parents 8f1eae1450b2
children e89715fe6a6e
comparison
equal deleted inserted replaced
699:8f1eae1450b2 700:de7f665e1d00
207 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 207 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
208 % type is a string specifying the type of boundary condition if there are several. 208 % type is a string specifying the type of boundary condition if there are several.
209 % data is a function returning the data that should be applied at the boundary. 209 % data is a function returning the data that should be applied at the boundary.
210 % neighbour_scheme is an instance of Scheme that should be interfaced to. 210 % neighbour_scheme is an instance of Scheme that should be interfaced to.
211 % neighbour_boundary is a string specifying which boundary to interface to. 211 % neighbour_boundary is a string specifying which boundary to interface to.
212 function [closure, penalty] = boundary_condition(obj, boundary,~) 212 function [closure, penalty,closureHamiltonian,penaltyHamiltonian] = boundary_condition(obj, boundary,~)
213 [e, d_n, d_t, coeff_t, coeff_n s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g] = obj.get_boundary_ops(boundary); 213 [e, d_n, d_t, coeff_t, coeff_n s, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t,g] = obj.get_boundary_ops(boundary);
214 214
215 a_t = @(t) spdiag(coeff_t(t)); 215 a_t = @(t) spdiag(coeff_t(t));
216 a_n = @(t) spdiag(coeff_n(t)); 216 a_n = @(t) spdiag(coeff_n(t));
217 217
224 penalty_parameter_1 = @(t) 1*1i*halfnorm_inv_n*halfnorm_inv_t*F(t)*e'*halfnorm_t*e; 224 penalty_parameter_1 = @(t) 1*1i*halfnorm_inv_n*halfnorm_inv_t*F(t)*e'*halfnorm_t*e;
225 penalty_parameter_2 = @(t) halfnorm_inv_n*e*tau2(t); 225 penalty_parameter_2 = @(t) halfnorm_inv_n*e*tau2(t);
226 226
227 closure = @(t) sqrt(obj.Ji)*(obj.c^2 * penalty_parameter_1(t)*e' + penalty_parameter_2(t)*e')*sqrt(obj.Ji); 227 closure = @(t) sqrt(obj.Ji)*(obj.c^2 * penalty_parameter_1(t)*e' + penalty_parameter_2(t)*e')*sqrt(obj.Ji);
228 penalty = @(t) -sqrt(obj.Ji)*(obj.c^2 * penalty_parameter_1(t)*e' + penalty_parameter_2(t)*e')*sqrt(obj.Ji); 228 penalty = @(t) -sqrt(obj.Ji)*(obj.c^2 * penalty_parameter_1(t)*e' + penalty_parameter_2(t)*e')*sqrt(obj.Ji);
229
230 closureHamiltonian = @(t) sqrt(obj.Ji)*(obj.c^2 * penalty_parameter_1(t)*e')*sqrt(obj.Ji);
231 penaltyHamiltonian = @(t) -sqrt(obj.Ji)*(obj.c^2 * penalty_parameter_1(t)*e')*sqrt(obj.Ji);
229 232
230 end 233 end
231 234
232 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 235 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
233 % u denotes the solution in the own domain 236 % u denotes the solution in the own domain