Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst2dCurve.m @ 1033:037f203b9bf5 feature/burgers1d

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with branch feature/advectioRV to utilize the +rv package
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 17 Jan 2019 10:44:12 +0100
parents 706d1c2b4199
children 78db023a7fe3
comparison
equal deleted inserted replaced
854:18162a0a5bb5 1033:037f203b9bf5
2 properties 2 properties
3 m % Number of points in each direction, possibly a vector 3 m % Number of points in each direction, possibly a vector
4 n % size of system 4 n % size of system
5 h % Grid spacing 5 h % Grid spacing
6 X,Y % Values of x and y for each grid point 6 X,Y % Values of x and y for each grid point
7 7
8 J, Ji % Jacobaian and inverse Jacobian 8 J, Ji % Jacobaian and inverse Jacobian
9 xi,eta 9 xi,eta
10 Xi,Eta 10 Xi,Eta
11 11
12 A,B 12 A,B
13 X_eta, Y_eta 13 X_eta, Y_eta
14 X_xi,Y_xi 14 X_xi,Y_xi
15 order % Order accuracy for the approximation 15 order % Order accuracy for the approximation
16 16
17 D % non-stabalized scheme operator 17 D % non-stabalized scheme operator
18 Ahat, Bhat, E 18 Ahat, Bhat, E
19 19
20 H % Discrete norm 20 H % Discrete norm
21 Hxii,Hetai % Kroneckerd norms in xi and eta. 21 Hxii,Hetai % Kroneckerd norms in xi and eta.
22 I_xi,I_eta, I_N, onesN 22 I_xi,I_eta, I_N, onesN
23 e_w, e_e, e_s, e_n 23 e_w, e_e, e_s, e_n
24 index_w, index_e,index_s,index_n 24 index_w, index_e,index_s,index_n
25 params % Parameters for the coeficient matrice 25 params % Parameters for the coeficient matrice
26 end 26 end
27 27
28 28
29 methods 29 methods
30 % Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Eu 30 % Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Eu
31 function obj = Hypsyst2dCurve(m, order, A, B, E, params, ti) 31 function obj = Hypsyst2dCurve(m, order, A, B, E, params, ti)
32 default_arg('E', []) 32 default_arg('E', [])
33 xilim = {0 1}; 33 xilim = {0 1};
34 etalim = {0 1}; 34 etalim = {0 1};
35 35
36 if length(m) == 1 36 if length(m) == 1
37 m = [m m]; 37 m = [m m];
38 end 38 end
39 obj.params = params; 39 obj.params = params;
40 obj.A=A; 40 obj.A=A;
41 obj.B=B; 41 obj.B=B;
42 42
43 obj.Ahat=@(params,x,y,x_eta,y_eta)(A(params,x,y).*y_eta-B(params,x,y).*x_eta); 43 obj.Ahat=@(params,x,y,x_eta,y_eta)(A(params,x,y).*y_eta-B(params,x,y).*x_eta);
44 obj.Bhat=@(params,x,y,x_xi,y_xi)(B(params,x,y).*x_xi-A(params,x,y).*y_xi); 44 obj.Bhat=@(params,x,y,x_xi,y_xi)(B(params,x,y).*x_xi-A(params,x,y).*y_xi);
45 obj.E=@(params,x,y,~,~)E(params,x,y); 45 obj.E=@(params,x,y,~,~)E(params,x,y);
46 46
47 m_xi = m(1); 47 m_xi = m(1);
48 m_eta = m(2); 48 m_eta = m(2);
49 m_tot=m_xi*m_eta; 49 m_tot=m_xi*m_eta;
50 50
51 ops_xi = sbp.D2Standard(m_xi,xilim,order); 51 ops_xi = sbp.D2Standard(m_xi,xilim,order);
52 ops_eta = sbp.D2Standard(m_eta,etalim,order); 52 ops_eta = sbp.D2Standard(m_eta,etalim,order);
53 53
54 obj.xi = ops_xi.x; 54 obj.xi = ops_xi.x;
55 obj.eta = ops_eta.x; 55 obj.eta = ops_eta.x;
56 56
57 obj.Xi = kr(obj.xi,ones(m_eta,1)); 57 obj.Xi = kr(obj.xi,ones(m_eta,1));
58 obj.Eta = kr(ones(m_xi,1),obj.eta); 58 obj.Eta = kr(ones(m_xi,1),obj.eta);
59 59
60 obj.n = length(A(obj.params,0,0)); 60 obj.n = length(A(obj.params,0,0));
61 obj.onesN=ones(obj.n); 61 obj.onesN=ones(obj.n);
62 62
63 obj.index_w=1:m_eta; 63 obj.index_w=1:m_eta;
64 obj.index_e=(m_tot-m_e 64 obj.index_e=(m_tot-m_e
65 65
66 metric_termsta+1):m_tot; 66 metric_termsta+1):m_tot;
67 obj.index_s=1:m_eta:(m_tot-m_eta+1); 67 obj.index_s=1:m_eta:(m_tot-m_eta+1);
68 obj.index_n=(m_eta):m_eta:m_tot; 68 obj.index_n=(m_eta):m_eta:m_tot;
69 69
70 I_n = eye(obj.n); 70 I_n = eye(obj.n);
71 I_xi = speye(m_xi); 71 I_xi = speye(m_xi);
72 obj.I_xi = I_xi; 72 obj.I_xi = I_xi;
73 I_eta = speye(m_eta); 73 I_eta = speye(m_eta);
74 obj.I_eta = I_eta; 74 obj.I_eta = I_eta;
75 75
76 D1_xi = kr(I_n, ops_xi.D1, I_eta); 76 D1_xi = kr(I_n, ops_xi.D1, I_eta);
77 obj.Hxii = kr(I_n, ops_xi.HI, I_eta); 77 obj.Hxii = kr(I_n, ops_xi.HI, I_eta);
78 D1_eta = kr(I_n, I_xi, ops_eta.D1); 78 D1_eta = kr(I_n, I_xi, ops_eta.D1);
79 obj.Hetai = kr(I_n, I_xi, ops_eta.HI); 79 obj.Hetai = kr(I_n, I_xi, ops_eta.HI);
80 80
81 obj.e_w = kr(I_n, ops_xi.e_l, I_eta); 81 obj.e_w = kr(I_n, ops_xi.e_l, I_eta);
82 obj.e_e = kr(I_n, ops_xi.e_r, I_eta); 82 obj.e_e = kr(I_n, ops_xi.e_r, I_eta);
83 obj.e_s = kr(I_n, I_xi, ops_eta.e_l); 83 obj.e_s = kr(I_n, I_xi, ops_eta.e_l);
84 obj.e_n = kr(I_n, I_xi, 84 obj.e_n = kr(I_n, I_xi,
85 85
86 metric_termsops_eta.e_r); 86 metric_termsops_eta.e_r);
87 87
88 [X,Y] = ti.map(obj.xi,obj.eta); 88 [X,Y] = ti.map(obj.xi,obj.eta);
89 89
90 [x_xi,x_eta] = gridDerivatives(X,ops_xi.D1,ops_eta.D1); 90 [x_xi,x_eta] = gridDerivatives(X,ops_xi.D1,ops_eta.D1);
91 [y_xi,y_eta] = gridDerivatives(Y,ops_xi.D1, ops_eta.D1); 91 [y_xi,y_eta] = gridDerivatives(Y,ops_xi.D1, ops_eta.D1);
92 92
93 obj.X = reshape(X,m_tot,1); 93 obj.X = reshape(X,m_tot,1);
94 obj.Y = reshape(Y,m_tot,1); 94 obj.Y = reshape(Y,m_tot,1);
95 obj.X_xi = reshape(x_xi,m_tot,1); 95 obj.X_xi = reshape(x_xi,m_tot,1);
96 obj.Y_xi = reshape(y_xi,m_tot,1); 96 obj.Y_xi = reshape(y_xi,m_tot,1);
97 obj.X_eta = reshape(x_eta,m_tot,1); 97 obj.X_eta = reshape(x_eta,m_tot,1);
98 obj.Y_eta = reshape(y_eta,m_tot,1); 98 obj.Y_eta = reshape(y_eta,m_tot,1);
99 99
100 Ahat_evaluated = obj.evaluateCoefficientMatrix(obj.Ahat, obj.X, obj.Y,obj.X_eta,obj.Y_eta); 100 Ahat_evaluated = obj.evaluateCoefficientMatrix(obj.Ahat, obj.X, obj.Y,obj.X_eta,obj.Y_eta);
101 Bhat_evaluated = obj.evaluateCoefficientMatrix(obj.Bhat, obj.X, obj.Y,obj.X_xi,obj.Y_xi); 101 Bhat_evaluated = obj.evaluateCoefficientMatrix(obj.Bhat, obj.X, obj.Y,obj.X_xi,obj.Y_xi);
102 E_evaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,[],[]); 102 E_evaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,[],[]);
103 103
104 obj.m = m; 104 obj.m = m;
105 obj.h = [ops_xi.h ops_eta.h]; 105 obj.h = [ops_xi.h ops_eta.h];
106 obj.order = order; 106 obj.order = order;
107 obj.J = obj.X_xi.*obj.Y_eta - obj.X_eta.*obj.Y_xi; 107 obj.J = obj.X_xi.*obj.Y_eta - obj.X_eta.*obj.Y_xi;
108 obj.Ji = kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot)); 108 obj.Ji = kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot));
109 109
110 obj.D = obj.Ji*(-Ahat_evaluated*D1_xi-Bhat_evaluated*D1_eta)-E_evaluated; 110 obj.D = obj.Ji*(-Ahat_evaluated*D1_xi-Bhat_evaluated*D1_eta)-E_evaluated;
111 111
112 end 112 end
113 113
114 % Closure functions return the opertors applied to the own doamin to close the boundary 114 % Closure functions return the opertors applied to the own doamin to close the boundary
115 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 115 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
116 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w',General boundary conditions'n','s'. 116 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w',General boundary conditions'n','s'.
117 % type is a string specifying the type of boundary condition if there are several. 117 % type is a string specifying the type of boundary condition if there are several.
118 % data is a function returning the data that should be applied at the boundary. 118 % data is a function returning the data that should be applied at the boundary.
125 [closure,penalty] = boundary_condition_general(obj,boundary,L); 125 [closure,penalty] = boundary_condition_general(obj,boundary,L);
126 otherwise 126 otherwise
127 error('No such boundary condition') 127 error('No such boundary condition')
128 end 128 end
129 end 129 end
130 130
131 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundaryGeneral boundary conditions) 131 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type)
132 error('An interface function does not exist yet'); 132 error('Not implemented');
133 end 133 end
134 134
135 function N = size(obj) 135 function N = size(obj)
136 N = obj.m; 136 N = obj.m;
137 end 137 end
138 138
139 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y,x_,y_) 139 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y,x_,y_)
140 params = obj.params; 140 params = obj.params;
141 141
142 if isa(mat,'function_handle') 142 if isa(mat,'function_handle')
143 [rows,cols] = size(mat(params,0,0,0,0)); 143 [rows,cols] = size(mat(params,0,0,0,0));
144 x_ = kr(obj.onesN,x_); 144 x_ = kr(obj.onesN,x_);
145 y_ = kr(obj.onesN,y_); 145 y_ = kr(obj.onesN,y_);
146 matVec = mat(params,X',Y',x_',y_'); 146 matVec = mat(params,X',Y',x_',y_');
150 matVec = mat; 150 matVec = mat;
151 [rows,cols] = size(matVec); 151 [rows,cols] = size(matVec);
152 side = max(length(X),length(Y)); 152 side = max(length(X),length(Y));
153 cols = cols/side; 153 cols = cols/side;
154 end 154 end
155 155
156 ret = cell(rows,cols); 156 ret = cell(rows,cols);
157 for ii = 1:rows 157 for ii = 1:rows
158 for jj = 1:cols 158 for jj = 1:cols
159 ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side)); 159 ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side));
160 end 160 end
161 end 161 end
162 ret = cell2mat(ret); 162 ret = cell2mat(ret);
163 end 163 end
164 164
165 %Characteristic boundary conditions 165 %Characteristic boundary conditions
166 function [closure, penalty] = boundary_condition_char(obj,boundary) 166 function [closure, penalty] = boundary_condition_char(obj,boundary)
167 params = obj.params; 167 params = obj.params;
168 X = obj.X; 168 X = obj.X;
169 Y = obj.Y; 169 Y = obj.Y;
170 xi = obj.xi; 170 xi = obj.xi;
171 eta = obj.eta; 171 eta = obj.eta;
172 172
173 switch boundary 173 switch boundary
174 case {'w','W','west'} 174 case {'w','W','west'}
175 e_ = obj.e_w; 175 e_ = obj.e_w;
176 mat = obj.Ahat; 176 mat = obj.Ahat;
177 boundPos = 'l'; 177 boundPos = 'l';
198 boundPos = 'r'; 198 boundPos = 'r';
199 Hi = obj.Hetai; 199 Hi = obj.Hetai;
200 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n)); 200 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n));
201 side = max(length(xi)); 201 side = max(length(xi));
202 end 202 end
203 203
204 pos = signVec(1); 204 pos = signVec(1);
205 zeroval = signVec(2); 205 zeroval = signVec(2);
206 neg = signVec(3); 206 neg = signVec(3);
207 207
208 switch boundPos 208 switch boundPos
209 case {'l'} 209 case {'l'}
210 tau = sparse(obj.n*side,pos); 210 tau = sparse(obj.n*side,pos);
211 Vi_plus = Vi(1:pos,:); 211 Vi_plus = Vi(1:pos,:);
212 tau(1:pos,:) = -abs(D(1:pos,1:pos)); 212 tau(1:pos,:) = -abs(D(1:pos,1:pos));
216 tau = sparse(obj.n*side,neg); 216 tau = sparse(obj.n*side,neg);
217 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); 217 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
218 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); 218 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
219 closure = Hi*e_*V*tau*Vi_minus*e_'; 219 closure = Hi*e_*V*tau*Vi_minus*e_';
220 penalty = -Hi*e_*V*tau*Vi_minus; 220 penalty = -Hi*e_*V*tau*Vi_minus;
221 end 221 end
222 end 222 end
223 223
224 224
225 % General boundary condition in the form Lu=g(x) 225 % General boundary condition in the form Lu=g(x)
226 function [closure,penalty] = boundary_condition_general(obj,boundary,L) 226 function [closure,penalty] = boundary_condition_general(obj,boundary,L)
227 params = obj.params; 227 params = obj.params;
228 X = obj.X; 228 X = obj.X;
229 Y = obj.Y; 229 Y = obj.Y;
230 xi = obj.xi; 230 xi = obj.xi;
231 eta = obj.eta; 231 eta = obj.eta;
232 232
233 switch boundary 233 switch boundary
234 case {'w','W','west'} 234 case {'w','W','west'}
235 e_ = obj.e_w; 235 e_ = obj.e_w;
236 mat = obj.Ahat; 236 mat = obj.Ahat;
237 boundPos = 'l'; 237 boundPos = 'l';
238 Hi = obj.Hxii; 238 Hi = obj.Hxii;
239 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w)); 239 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w));
240 240
241 Ji_vec = diag(obj.Ji); 241 Ji_vec = diag(obj.Ji);
242 Ji = diag(Ji_vec(obj.index_w)); 242 Ji = diag(Ji_vec(obj.index_w));
243 xi_x = Ji*obj.Y_eta(obj.index_w); 243 xi_x = Ji*obj.Y_eta(obj.index_w);
244 xi_y = -Ji*obj.X_eta(obj.index_w); 244 xi_y = -Ji*obj.X_eta(obj.index_w);
245 L = obj.evaluateCoefficientMatrix(L,xi_x,xi_y,[],[]); 245 L = obj.evaluateCoefficientMatrix(L,xi_x,xi_y,[],[]);
248 e_ = obj.e_e; 248 e_ = obj.e_e;
249 mat = obj.Ahat; 249 mat = obj.Ahat;
250 boundPos = 'r'; 250 boundPos = 'r';
251 Hi = obj.Hxii; 251 Hi = obj.Hxii;
252 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e)); 252 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e));
253 253
254 Ji_vec = diag(obj.Ji); 254 Ji_vec = diag(obj.Ji);
255 Ji = diag(Ji_vec(obj.index_e)); 255 Ji = diag(Ji_vec(obj.index_e));
256 xi_x = Ji*obj.Y_eta(obj.index_e); 256 xi_x = Ji*obj.Y_eta(obj.index_e);
257 xi_y = -Ji*obj.X_eta(obj.index_e); 257 xi_y = -Ji*obj.X_eta(obj.index_e);
258 L = obj.evaluateCoefficientMatrix(L,-xi_x,-xi_y,[],[]); 258 L = obj.evaluateCoefficientMatrix(L,-xi_x,-xi_y,[],[]);
261 e_ = obj.e_s; 261 e_ = obj.e_s;
262 mat = obj.Bhat; 262 mat = obj.Bhat;
263 boundPos = 'l'; 263 boundPos = 'l';
264 Hi = obj.Hetai; 264 Hi = obj.Hetai;
265 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s)); 265 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s));
266 266
267 Ji_vec = diag(obj.Ji); 267 Ji_vec = diag(obj.Ji);
268 Ji = diag(Ji_vec(obj.index_s)); 268 Ji = diag(Ji_vec(obj.index_s));
269 eta_x = Ji*obj.Y_xi(obj.index_s); 269 eta_x = Ji*obj.Y_xi(obj.index_s);
270 eta_y = -Ji*obj.X_xi(obj.index_s); 270 eta_y = -Ji*obj.X_xi(obj.index_s);
271 L = obj.evaluateCoefficientMatrix(L,eta_x,eta_y,[],[]); 271 L = obj.evaluateCoefficientMatrix(L,eta_x,eta_y,[],[]);
274 e_ = obj.e_n; 274 e_ = obj.e_n;
275 mat = obj.Bhat; 275 mat = obj.Bhat;
276 boundPos = 'r'; 276 boundPos = 'r';
277 Hi = obj.Hetai; 277 Hi = obj.Hetai;
278 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n)); 278 [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n));
279 279
280 Ji_vec = diag(obj.Ji); 280 Ji_vec = diag(obj.Ji);
281 Ji = diag(Ji_vec(obj.index_n)); 281 Ji = diag(Ji_vec(obj.index_n));
282 eta_x = Ji*obj.Y_xi(obj.index_n); 282 eta_x = Ji*obj.Y_xi(obj.index_n);
283 eta_y = -Ji*obj.X_xi(obj.index_n); 283 eta_y = -Ji*obj.X_xi(obj.index_n);
284 L = obj.evaluateCoefficientMatrix(L,-eta_x,-eta_y,[],[]); 284 L = obj.evaluateCoefficientMatrix(L,-eta_x,-eta_y,[],[]);
285 side = max(length(xi)); 285 side = max(length(xi));
286 end 286 end
287 287
288 pos = signVec(1); 288 pos = signVec(1);
289 zeroval = signVec(2); 289 zeroval = signVec(2);
290 neg = signVec(3); 290 neg = signVec(3);
291 291
292 switch boundPos 292 switch boundPos
293 case {'l'} 293 case {'l'}
294 tau = sparse(obj.n*side,pos); 294 tau = sparse(obj.n*side,pos);
295 Vi_plus = Vi(1:pos,:); 295 Vi_plus = Vi(1:pos,:);
296 Vi_minus = Vi(pos+1:obj.n*side,:); 296 Vi_minus = Vi(pos+1:obj.n*side,:);
297 V_plus = V(:,1:pos); 297 V_plus = V(:,1:pos);
298 V_minus = V(:,(pos)+1:obj.n*side); 298 V_minus = V(:,(pos)+1:obj.n*side);
299 299
300 tau(1:pos,:) = -abs(D(1:pos,1:pos)); 300 tau(1:pos,:) = -abs(D(1:pos,1:pos));
301 R = -inv(L*V_plus)*(L*V_minus); 301 R = -inv(L*V_plus)*(L*V_minus);
302 closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; 302 closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
303 penalty = -Hi*e_*V*tau*inv(L*V_plus)*L; 303 penalty = -Hi*e_*V*tau*inv(L*V_plus)*L;
304 case {'r'} 304 case {'r'}
305 tau = sparse(obj.n*side,neg); 305 tau = sparse(obj.n*side,neg);
306 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); 306 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
307 Vi_plus = Vi(1:pos,:); 307 Vi_plus = Vi(1:pos,:);
308 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); 308 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
309 309
310 V_plus = V(:,1:pos); 310 V_plus = V(:,1:pos);
311 V_minus = V(:,(pos+zeroval)+1:obj.n*side); 311 V_minus = V(:,(pos+zeroval)+1:obj.n*side);
312 R = -inv(L*V_minus)*(L*V_plus); 312 R = -inv(L*V_minus)*(L*V_plus);
313 closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_'; 313 closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
314 penalty = -Hi*e_*V*tau*inv(L*V_minus)*L; 314 penalty = -Hi*e_*V*tau*inv(L*V_minus)*L;
315 end 315 end
316 end 316 end
317 317
318 % Function that diagonalizes a symbolic matrix A as A=V*D*Vi 318 % Function that diagonalizes a symbolic matrix A as A=V*D*Vi
319 % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign 319 % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign
320 % [d+ ] 320 % [d+ ]
321 % D = [ d0 ] 321 % D = [ d0 ]
322 % [ d-] 322 % [ d-]
327 if(sum(abs(x_)) ~= 0) 327 if(sum(abs(x_)) ~= 0)
328 syms xs_ 328 syms xs_
329 else 329 else
330 xs_ = 0; 330 xs_ = 0;
331 end 331 end
332 332
333 if(sum(abs(y_))~= 0) 333 if(sum(abs(y_))~= 0)
334 syms ys_; 334 syms ys_;
335 else 335 else
336 ys_ = 0; 336 ys_ = 0;
337 end 337 end
338 338
339 [V, D] = eig(mat(params,xs,ys,xs_,ys_)); 339 [V, D] = eig(mat(params,xs,ys,xs_,ys_));
340 Vi = inv(V); 340 Vi = inv(V);
341 syms xs ys xs_ ys_ 341 syms xs ys xs_ ys_
342 342
343 xs = x; 343 xs = x;
344 ys = y; 344 ys = y;
345 xs_ = x_; 345 xs_ = x_;
346 ys_ = y_; 346 ys_ = y_;
347 347
348 side = max(length(x),length(y)); 348 side = max(length(x),length(y));
349 Dret = zeros(obj.n,side*obj.n); 349 Dret = zeros(obj.n,side*obj.n);
350 Vret = zeros(obj.n,side*obj.n); 350 Vret = zeros(obj.n,side*obj.n);
351 Viret = zeros(obj.n,side*obj.n); 351 Viret = zeros(obj.n,side*obj.n);
352 for ii = 1:obj.n 352 for ii = 1:obj.n
354 Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii)); 354 Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii));
355 Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii)); 355 Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii));
356 Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii)); 356 Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii));
357 end 357 end
358 end 358 end
359 359
360 D = sparse(Dret); 360 D = sparse(Dret);
361 V = sparse(Vret); 361 V = sparse(Vret);
362 Vi = sparse(Viret); 362 Vi = sparse(Viret);
363 V = obj.evaluateCoefficientMatrix(V,x,y,x_,y_); 363 V = obj.evaluateCoefficientMatrix(V,x,y,x_,y_);
364 D = obj.evaluateCoefficientMatrix(D,x,y,x_,y_); 364 D = obj.evaluateCoefficientMatrix(D,x,y,x_,y_);
365 Vi = obj.evaluateCoefficientMatrix(Vi,x,y,x_,y_); 365 Vi = obj.evaluateCoefficientMatrix(Vi,x,y,x_,y_);
366 DD = diag(D); 366 DD = diag(D);
367 367
368 poseig = (DD>0); 368 poseig = (DD>0);
369 zeroeig = (DD==0); 369 zeroeig = (DD==0);
370 negeig = (DD<0); 370 negeig = (DD<0);
371 371
372 D = diag([DD(poseig); DD(zeroeig); DD(negeig)]); 372 D = diag([DD(poseig); DD(zeroeig); DD(negeig)]);
373 V = [V(:,poseig) V(:,zeroeig) V(:,negeig)]; 373 V = [V(:,poseig) V(:,zeroeig) V(:,negeig)];
374 Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)]; 374 Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
375 signVec = [sum(poseig),sum(zeroeig),sum(negeig)]; 375 signVec = [sum(poseig),sum(zeroeig),sum(negeig)];
376 end 376 end