Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger1dCurve.m @ 508:508b7493be94 feature/quantumTriangles

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
moved all multiplications and divisions bu J and Ji to the siscretization to easier find the error
author Ylva Rydin <ylva.rydin@telia.com>
date Fri, 02 Jun 2017 11:33:04 +0200
parents 1341e8580d64
children 57f3493f851b
comparison
equal deleted inserted replaced
507:1341e8580d64 508:508b7493be94
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 = (-0.5*(obj.D1*obj.a + obj.a*obj.D1) + 1i*obj.D2(diag(obj.Ji)) + 1i*obj.V_mat);
89 90 % D = (-0.5*(obj.D1*obj.a -obj.a_xi+ obj.a*obj.D1) + 1i*obj.D2(diag(obj.Ji)) + 1i*obj.V_mat);
91 % 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 92 end
91 93
92 94
93 function [] = variable_update(obj,t) 95 function [] = variable_update(obj,t)
94 if (t == obj.t_up) 96 if (t == obj.t_up)
100 ddt_x_l = obj.ddt_x_l(t); 102 ddt_x_l = obj.ddt_x_l(t);
101 obj.x = obj.xi*(x_r -x_l) + x_l; 103 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))); 104 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 105
104 obj.Ji = sparse(diag(1./(x_r - x_l + 0*obj.x))); 106 obj.Ji = sparse(diag(1./(x_r - x_l + 0*obj.x)));
107 obj.J = sparse(x_r -x_l);
105 obj.a_xi = sparse(diag(-1*(ddt_x_r - ddt_x_l + 0*obj.x))); 108 obj.a_xi = sparse(diag(-1*(ddt_x_r - ddt_x_l + 0*obj.x)));
106 obj.t_up = t; 109 obj.t_up = t;
107 end 110 end
108 end 111 end
109 112
114 [e,d,s,p] = obj.get_boundary_ops(boundary); 117 [e,d,s,p] = obj.get_boundary_ops(boundary);
115 118
116 switch type 119 switch type
117 % Dirichlet boundary condition 120 % Dirichlet boundary condition
118 case {'D','d','dirichlet'} 121 case {'D','d','dirichlet'}
119 tau1 = @(t) s * 1i*obj.Ji(p,p)^2*d; 122 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; 123 tau2 = @(t) (1*s*obj.a(p,p))/2*e;
121 closure = @(t) obj.Hi*tau1(t)*e' + obj.Hi*tau2(obj.a)*e'; 124 closure = @(t)(obj.Hi*tau1(t)*e' + obj.Hi*tau2(obj.a)*e');
122 125
123 switch class(data) 126 switch class(data)
124 case 'double' 127 case 'double'
125 penalty = @(t) -obj.Ji*(obj.Hi*tau1*data+obj.Hi*tau2(obj.a)*data); 128 penalty = @(t) -(obj.Hi*tau1*data+obj.Hi*tau2(obj.a)*data);
126 % case 'function_handle' 129 % case 'function_handle'
127 % penalty = @(t)-obj.Hi*tau*data(t); 130 % penalty = @(t)-obj.Hi*tau*data(t);
128 otherwise 131 otherwise
129 error('Wierd data argument!') 132 error('Weird data argument!')
130 end 133 end
131 134
132 % Unknown, boundary condition 135 % Unknown, boundary condition
133 otherwise 136 otherwise
134 error('No such boundary condition: type = %s',type); 137 error('No such boundary condition: type = %s',type);
141 [e_u,d_u,s_u,p_u] = obj.get_boundary_ops(boundary); 144 [e_u,d_u,s_u,p_u] = obj.get_boundary_ops(boundary);
142 [e_v,d_v,s_v,p_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 145 [e_v,d_v,s_v,p_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
143 146
144 a1 = s_u* 1/2 * 1i ; 147 a1 = s_u* 1/2 * 1i ;
145 b1 = -s_u* 1/2 * 1i; 148 b1 = -s_u* 1/2 * 1i;
146 gamma = @(a) -obj.Ji*s_u*a(p_u,p_u)/2*e_u; 149 gamma = @(a) -s_u*a(p_u,p_u)/2*e_u;
147 150
148 tau = @(t) a1*obj.Ji(p_u,p_u)^2*d_u; 151 tau = @(t) a1*obj.Ji(p_u,p_u)^2*d_u;
149 sig = b1*e_u; 152 sig = b1*e_u;
150 153
151 closure = @(t) obj.Hi * (tau(t)*e_u' + sig*obj.Ji(p_u,p_u)^2*d_u' + gamma(obj.a)*e_u'); 154 closure = @(t) obj.Hi * (tau(t)*e_u' + sig*obj.Ji(p_u,p_u)^2*d_u' + gamma(obj.a)*e_u');