Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger1dCurve.m @ 496:437fad4a47e1 feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add 1d interface for shrodinger in 1d
author Ylva Rydin <ylva.rydin@telia.com>
date Fri, 24 Feb 2017 13:58:18 +0100
parents b91d23271481
children 324c927d8b1d
comparison
equal deleted inserted replaced
495:b91d23271481 496:437fad4a47e1
8 8
9 D % non-stabalized scheme operator 9 D % non-stabalized scheme operator
10 H % Discrete norm 10 H % Discrete norm
11 M % Derivative norm 11 M % Derivative norm
12 alpha 12 alpha
13 x_r
14 x_l
15 ddt_x_r
16 ddt_x_l
17 a
18 a_xi
19 Ji
20 t_up
21 x
13 22
14 V_mat 23 V_mat
15 D1 24 D1
16 D2 25 D2
17 Hi 26 Hi
22 gamm 31 gamm
23 end 32 end
24 33
25 methods 34 methods
26 % Solving SE in the form u_t = i*u_xx +i*V on deforming 1D domain; 35 % Solving SE in the form u_t = i*u_xx +i*V on deforming 1D domain;
27 function obj = Schrodinger1dCurve(g,m,order,V,constJi) 36 function obj = Schrodinger1dCurve(g,order,boundaries,V,constJi)
28 default_arg('V',0); 37 default_arg('V',0);
29 default_arg('constJi',false) 38 default_arg('constJi',false)
30 xilim={0 1}; 39 xilim={0 1};
40 m = N(g);
31 if constJi 41 if constJi
32 ops = sbp.D2Standard(m,xilim,order); 42 ops = sbp.D2Standard(m,xilim,order);
33 else 43 else
34 ops = sbp.D4Variable(m,xilim,order); 44 ops = sbp.D4Variable(m,xilim,order);
35 end 45 end
46
47 obj.x_l = boundaries{1};
48 obj.x_r = boundaries{2};
49 obj.ddt_x_l = boundaries{3};
50 obj.ddt_x_r = boundaries{4};
36 51
37 obj.xi=ops.x; 52 obj.xi=ops.x;
38 obj.h=ops.h; 53 obj.h=ops.h;
39 obj.D2 = ops.D2; 54 obj.D2 = ops.D2;
40 obj.D1 = ops.D1; 55 obj.D1 = ops.D1;
45 obj.e_r = ops.e_r; 60 obj.e_r = ops.e_r;
46 obj.d1_l = ops.d1_l; 61 obj.d1_l = ops.d1_l;
47 obj.d1_r = ops.d1_r; 62 obj.d1_r = ops.d1_r;
48 obj.grid = g; 63 obj.grid = g;
49 64
50
51 if isa(V,'function_handle') 65 if isa(V,'function_handle')
52 V_vec = V(obj.x); 66 V_vec = V(obj.x);
53 else 67 else
54 V_vec = obj.xi*0 + V; 68 V_vec = obj.xi*0 + V;
55 end 69 end
56 70
57 obj.V_mat = spdiags(V_vec,0,m,m); 71 obj.V_mat = spdiags(V_vec,0,m,m);
58 obj.D = @(a,a_xi,Ji) obj.d_fun(a, a_xi, Ji, constJi); 72 obj.D = @(t) obj.d_fun(t);
59 obj.m = m; 73 obj.m = m;
60 obj.order = order; 74 obj.order = order;
61 end 75 end
62 76
63 77
67 % type is a string specifying the type of boundary condition if there are several. 81 % type is a string specifying the type of boundary condition if there are several.
68 % data is a function returning the data that should be applied at the boundary. 82 % data is a function returning the data that should be applied at the boundary.
69 % 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.
70 % neighbour_boundary is a string specifying which boundary to interface to. 84 % neighbour_boundary is a string specifying which boundary to interface to.
71 85
72 function [D] = d_fun(obj,a, a_xi , Ji , constJi) 86 function [D] = d_fun(obj,t)
73 if constJi 87 % obj.update_vairables(t); In driscretization?
74 D= -0.5*(obj.D1*a - a_xi + a*obj.D1) + 1i*Ji*obj.D2 + 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);
89
90 end
91
92
93 function [] = variable_update(obj,t)
94 if (t == obj.t_up)
95 return
75 else 96 else
76 D= -0.5*(obj.D1*a - a_xi + a*obj.D1) + 1i*obj.D2(diag(Ji)) + 1i*obj.V_mat; 97 x_r = obj.x_r(t);
98 x_l = obj.x_l(t);
99 ddt_x_r = obj. ddt_x_r(t);
100 ddt_x_l = obj.ddt_x_l(t);
101 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
104 obj.Ji = sparse(diag(1./(x_r - x_l + 0*obj.x)));
105 obj.a_xi = sparse(diag(-1*(ddt_x_r - ddt_x_l + 0*obj.x)));
106 obj.t_up = t;
77 end 107 end
78 end 108 end
79 109
80 function [closure, penalty] = boundary_condition(obj,boundary,type,data) 110 function [closure, penalty] = boundary_condition(obj,boundary,type,data)
81 default_arg('type','dirichlet'); 111 default_arg('type','dirichlet');
85 115
86 switch type 116 switch type
87 % Dirichlet boundary condition 117 % Dirichlet boundary condition
88 case {'D','d','dirichlet'} 118 case {'D','d','dirichlet'}
89 tau1 = s * 1i*d; 119 tau1 = s * 1i*d;
90 tau2 = @(a) (-1*s*a(p,p) - abs(a(p,p)))/4*e; 120 tau2 = @(t) obj.Ji*(-1*s*obj.a(p,p) - abs(obj.a(p,p)))/4*e;
91 closure = @(a) obj.Hi*tau1*e' + obj.Hi*tau2(a)*e'; 121 closure = @(t) obj.Ji*obj.Hi*tau1*e' + obj.Hi*tau2(obj.a)*e';
92 122
93 switch class(data) 123 switch class(data)
94 case 'double' 124 case 'double'
95 penalty = @(a) -(obj.Hi*tau1*data+obj.Hi*tau2(a)*data); 125 penalty = @(t) -obj.Ji*(obj.Hi*tau1*data+obj.Hi*tau2(obj.a)*data);
96 % case 'function_handle' 126 % case 'function_handle'
97 % penalty = @(t)-obj.Hi*tau*data(t); 127 % penalty = @(t)-obj.Hi*tau*data(t);
98 otherwise 128 otherwise
99 error('Wierd data argument!') 129 error('Wierd data argument!')
100 end 130 end
106 end 136 end
107 137
108 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 138 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
109 % u denotes the solution in the own domain 139 % u denotes the solution in the own domain
110 % v denotes the solution in the neighbour domain 140 % v denotes the solution in the neighbour domain
111 [e_u,d_u,s_u,p_u] = obj.get_boundary_ops(boundary); 141 [e_u,d_u,s_u,p_u] = obj.get_boundary_ops(boundary);
112 [e_v,d_v,s_v,p_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 142 [e_v,d_v,s_v,p_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
113 143
114 a1 = -s_u* 1/2 * 1i ; 144 a1 = -s_u* 1/2 * 1i ;
115 b1 = a1'; 145 b1 = a1';
116 gamma = -1*s*a(p_u,p_u)/2; 146 gamma = @(a) -s_u*a(p_u,p_u)/2*e_u;
117 147
118 tau = b1*d_u; 148 tau = b1*d_u;
119 sig = -a1*e_u; 149 sig = -a1*e_u;
120 150
121 closure = @(a) obj.Hi * (tau*e_u' + sig*d_u' + gamma(a)); 151 closure = @(t) obj.Hi * (tau*e_u' + sig*d_u' + gamma(obj.a)*e_u');
122 penalty = @(a) obj.Hi * (-tau*e_v' - sig*d_v' - gamma(a)); 152 penalty = @(t) obj.Hi * (-tau*e_v' - sig*d_v' - gamma(obj.a)*e_v');
123 end 153 end
124 154
125 % Ruturns the boundary ops and sign for the boundary specified by the string boundary. 155 % Ruturns the boundary ops and sign for the boundary specified by the string boundary.
126 % The right boundary is considered the positive boundary 156 % The right boundary is considered the positive boundary
127 function [e,d,s,p] = get_boundary_ops(obj,boundary) 157 function [e,d,s,p] = get_boundary_ops(obj,boundary)