Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst3d.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 0652b34f9f27
comparison
equal deleted inserted replaced
854:18162a0a5bb5 1033:037f203b9bf5
5 h % Grid spacing 5 h % Grid spacing
6 x, y, z % Grid 6 x, y, z % Grid
7 X, Y, Z% Values of x and y for each grid point 7 X, Y, Z% Values of x and y for each grid point
8 Yx, Zx, Xy, Zy, Xz, Yz %Grid values for boundary surfaces 8 Yx, Zx, Xy, Zy, Xz, Yz %Grid values for boundary surfaces
9 order % Order accuracy for the approximation 9 order % Order accuracy for the approximation
10 10
11 D % non-stabalized scheme operator 11 D % non-stabalized scheme operator
12 A, B, C, E % Symbolic coefficient matrices 12 A, B, C, E % Symbolic coefficient matrices
13 Aevaluated,Bevaluated,Cevaluated, Eevaluated 13 Aevaluated,Bevaluated,Cevaluated, Eevaluated
14 14
15 H % Discrete norm 15 H % Discrete norm
16 Hx, Hy, Hz % Norms in the x, y and z directions 16 Hx, Hy, Hz % Norms in the x, y and z directions
17 Hxi,Hyi, Hzi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. 17 Hxi,Hyi, Hzi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
18 I_x,I_y, I_z, I_N 18 I_x,I_y, I_z, I_N
19 e_w, e_e, e_s, e_n, e_b, e_t 19 e_w, e_e, e_s, e_n, e_b, e_t
20 params % Parameters for the coeficient matrice 20 params % Parameters for the coeficient matrice
21 end 21 end
22 22
23 23
24 methods 24 methods
25 % Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Cu_z-Eu 25 % Solving Hyperbolic systems on the form u_t=-Au_x-Bu_y-Cu_z-Eu
26 function obj = Hypsyst3d(m, lim, order, A, B,C, E, params,operator) 26 function obj = Hypsyst3d(m, lim, order, A, B,C, E, params,operator)
27 default_arg('E', []) 27 default_arg('E', [])
28 xlim = lim{1}; 28 xlim = lim{1};
29 ylim = lim{2}; 29 ylim = lim{2};
30 zlim = lim{3}; 30 zlim = lim{3};
31 31
32 if length(m) == 1 32 if length(m) == 1
33 m = [m m m]; 33 m = [m m m];
34 end 34 end
35 35
36 obj.A = A; 36 obj.A = A;
37 obj.B = B; 37 obj.B = B;
38 obj.C = C; 38 obj.C = C;
39 obj.E = E; 39 obj.E = E;
40 m_x = m(1); 40 m_x = m(1);
41 m_y = m(2); 41 m_y = m(2);
42 m_z = m(3); 42 m_z = m(3);
43 obj.params = params; 43 obj.params = params;
44 44
45 switch operator 45 switch operator
46 case 'upwind' 46 case 'upwind'
47 ops_x = sbp.D1Upwind(m_x,xlim,order); 47 ops_x = sbp.D1Upwind(m_x,xlim,order);
48 ops_y = sbp.D1Upwind(m_y,ylim,order); 48 ops_y = sbp.D1Upwind(m_y,ylim,order);
49 ops_z = sbp.D1Upwind(m_z,zlim,order); 49 ops_z = sbp.D1Upwind(m_z,zlim,order);
50 otherwise 50 otherwise
51 ops_x = sbp.D2Standard(m_x,xlim,order); 51 ops_x = sbp.D2Standard(m_x,xlim,order);
52 ops_y = sbp.D2Standard(m_y,ylim,order); 52 ops_y = sbp.D2Standard(m_y,ylim,order);
53 ops_z = sbp.D2Standard(m_z,zlim,order); 53 ops_z = sbp.D2Standard(m_z,zlim,order);
54 end 54 end
55 55
56 obj.x = ops_x.x; 56 obj.x = ops_x.x;
57 obj.y = ops_y.x; 57 obj.y = ops_y.x;
58 obj.z = ops_z.x; 58 obj.z = ops_z.x;
59 59
60 obj.X = kr(obj.x,ones(m_y,1),ones(m_z,1)); 60 obj.X = kr(obj.x,ones(m_y,1),ones(m_z,1));
61 obj.Y = kr(ones(m_x,1),obj.y,ones(m_z,1)); 61 obj.Y = kr(ones(m_x,1),obj.y,ones(m_z,1));
62 obj.Z = kr(ones(m_x,1),ones(m_y,1),obj.z); 62 obj.Z = kr(ones(m_x,1),ones(m_y,1),obj.z);
63 63
64 obj.Yx = kr(obj.y,ones(m_z,1)); 64 obj.Yx = kr(obj.y,ones(m_z,1));
65 obj.Zx = kr(ones(m_y,1),obj.z); 65 obj.Zx = kr(ones(m_y,1),obj.z);
66 obj.Xy = kr(obj.x,ones(m_z,1)); 66 obj.Xy = kr(obj.x,ones(m_z,1));
67 obj.Zy = kr(ones(m_x,1),obj.z); 67 obj.Zy = kr(ones(m_x,1),obj.z);
68 obj.Xz = kr(obj.x,ones(m_y,1)); 68 obj.Xz = kr(obj.x,ones(m_y,1));
69 obj.Yz = kr(ones(m_z,1),obj.y); 69 obj.Yz = kr(ones(m_z,1),obj.y);
70 70
71 obj.Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y,obj.Z); 71 obj.Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y,obj.Z);
72 obj.Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y,obj.Z); 72 obj.Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y,obj.Z);
73 obj.Cevaluated = obj.evaluateCoefficientMatrix(C, obj.X, obj.Y,obj.Z); 73 obj.Cevaluated = obj.evaluateCoefficientMatrix(C, obj.X, obj.Y,obj.Z);
74 obj.Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y,obj.Z); 74 obj.Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y,obj.Z);
75 75
76 obj.n = length(A(obj.params,0,0,0)); 76 obj.n = length(A(obj.params,0,0,0));
77 77
78 I_n = speye(obj.n); 78 I_n = speye(obj.n);
79 I_x = speye(m_x); 79 I_x = speye(m_x);
80 obj.I_x = I_x; 80 obj.I_x = I_x;
81 I_y = speye(m_y); 81 I_y = speye(m_y);
82 obj.I_y = I_y; 82 obj.I_y = I_y;
83 I_z = speye(m_z); 83 I_z = speye(m_z);
84 obj.I_z = I_z; 84 obj.I_z = I_z;
85 I_N = kr(I_n,I_x,I_y,I_z); 85 I_N = kr(I_n,I_x,I_y,I_z);
86 86
87 obj.Hxi = kr(I_n, ops_x.HI, I_y,I_z); 87 obj.Hxi = kr(I_n, ops_x.HI, I_y,I_z);
88 obj.Hx = ops_x.H; 88 obj.Hx = ops_x.H;
89 obj.Hyi = kr(I_n, I_x, ops_y.HI,I_z); 89 obj.Hyi = kr(I_n, I_x, ops_y.HI,I_z);
90 obj.Hy = ops_y.H; 90 obj.Hy = ops_y.H;
91 obj.Hzi = kr(I_n, I_x,I_y, ops_z.HI); 91 obj.Hzi = kr(I_n, I_x,I_y, ops_z.HI);
92 obj.Hz = ops_z.H; 92 obj.Hz = ops_z.H;
93 93
94 obj.e_w = kr(I_n, ops_x.e_l, I_y,I_z); 94 obj.e_w = kr(I_n, ops_x.e_l, I_y,I_z);
95 obj.e_e = kr(I_n, ops_x.e_r, I_y,I_z); 95 obj.e_e = kr(I_n, ops_x.e_r, I_y,I_z);
96 obj.e_s = kr(I_n, I_x, ops_y.e_l,I_z); 96 obj.e_s = kr(I_n, I_x, ops_y.e_l,I_z);
97 obj.e_n = kr(I_n, I_x, ops_y.e_r,I_z); 97 obj.e_n = kr(I_n, I_x, ops_y.e_r,I_z);
98 obj.e_b = kr(I_n, I_x, I_y, ops_z.e_l); 98 obj.e_b = kr(I_n, I_x, I_y, ops_z.e_l);
99 obj.e_t = kr(I_n, I_x, I_y, ops_z.e_r); 99 obj.e_t = kr(I_n, I_x, I_y, ops_z.e_r);
100 100
101 obj.m = m; 101 obj.m = m;
102 obj.h = [ops_x.h ops_y.h ops_x.h]; 102 obj.h = [ops_x.h ops_y.h ops_x.h];
103 obj.order = order; 103 obj.order = order;
104 104
105 switch operator 105 switch operator
106 case 'upwind' 106 case 'upwind'
107 alphaA = max(abs(eig(A(params,obj.x(end),obj.y(end),obj.z(end))))); 107 alphaA = max(abs(eig(A(params,obj.x(end),obj.y(end),obj.z(end)))));
108 alphaB = max(abs(eig(B(params,obj.x(end),obj.y(end),obj.z(end))))); 108 alphaB = max(abs(eig(B(params,obj.x(end),obj.y(end),obj.z(end)))));
109 alphaC = max(abs(eig(C(params,obj.x(end),obj.y(end),obj.z(end))))); 109 alphaC = max(abs(eig(C(params,obj.x(end),obj.y(end),obj.z(end)))));
110 110
111 Ap = (obj.Aevaluated+alphaA*I_N)/2; 111 Ap = (obj.Aevaluated+alphaA*I_N)/2;
112 Am = (obj.Aevaluated-alphaA*I_N)/2; 112 Am = (obj.Aevaluated-alphaA*I_N)/2;
113 Dpx = kr(I_n, ops_x.Dp, I_y,I_z); 113 Dpx = kr(I_n, ops_x.Dp, I_y,I_z);
114 Dmx = kr(I_n, ops_x.Dm, I_y,I_z); 114 Dmx = kr(I_n, ops_x.Dm, I_y,I_z);
115 obj.D = -Am*Dpx; 115 obj.D = -Am*Dpx;
116 temp = Ap*Dmx; 116 temp = Ap*Dmx;
117 obj.D = obj.D-temp; 117 obj.D = obj.D-temp;
118 clear Ap Am Dpx Dmx 118 clear Ap Am Dpx Dmx
119 119
120 Bp = (obj.Bevaluated+alphaB*I_N)/2; 120 Bp = (obj.Bevaluated+alphaB*I_N)/2;
121 Bm = (obj.Bevaluated-alphaB*I_N)/2; 121 Bm = (obj.Bevaluated-alphaB*I_N)/2;
122 Dpy = kr(I_n, I_x, ops_y.Dp,I_z); 122 Dpy = kr(I_n, I_x, ops_y.Dp,I_z);
123 Dmy = kr(I_n, I_x, ops_y.Dm,I_z); 123 Dmy = kr(I_n, I_x, ops_y.Dm,I_z);
124 temp = Bm*Dpy; 124 temp = Bm*Dpy;
125 obj.D = obj.D-temp; 125 obj.D = obj.D-temp;
126 temp = Bp*Dmy; 126 temp = Bp*Dmy;
127 obj.D = obj.D-temp; 127 obj.D = obj.D-temp;
128 clear Bp Bm Dpy Dmy 128 clear Bp Bm Dpy Dmy
129 129
130 130
131 Cp = (obj.Cevaluated+alphaC*I_N)/2; 131 Cp = (obj.Cevaluated+alphaC*I_N)/2;
132 Cm = (obj.Cevaluated-alphaC*I_N)/2; 132 Cm = (obj.Cevaluated-alphaC*I_N)/2;
133 Dpz = kr(I_n, I_x, I_y,ops_z.Dp); 133 Dpz = kr(I_n, I_x, I_y,ops_z.Dp);
134 Dmz = kr(I_n, I_x, I_y,ops_z.Dm); 134 Dmz = kr(I_n, I_x, I_y,ops_z.Dm);
135 135
136 temp = Cm*Dpz; 136 temp = Cm*Dpz;
137 obj.D = obj.D-temp; 137 obj.D = obj.D-temp;
138 temp = Cp*Dmz; 138 temp = Cp*Dmz;
139 obj.D = obj.D-temp; 139 obj.D = obj.D-temp;
140 clear Cp Cm Dpz Dmz 140 clear Cp Cm Dpz Dmz
141 obj.D = obj.D-obj.Eevaluated; 141 obj.D = obj.D-obj.Eevaluated;
142 142
143 case 'standard' 143 case 'standard'
144 D1_x = kr(I_n, ops_x.D1, I_y,I_z); 144 D1_x = kr(I_n, ops_x.D1, I_y,I_z);
145 D1_y = kr(I_n, I_x, ops_y.D1,I_z); 145 D1_y = kr(I_n, I_x, ops_y.D1,I_z);
146 D1_z = kr(I_n, I_x, I_y,ops_z.D1); 146 D1_z = kr(I_n, I_x, I_y,ops_z.D1);
147 obj.D = -obj.Aevaluated*D1_x-obj.Bevaluated*D1_y-obj.Cevaluated*D1_z-obj.Eevaluated; 147 obj.D = -obj.Aevaluated*D1_x-obj.Bevaluated*D1_y-obj.Cevaluated*D1_z-obj.Eevaluated;
148 otherwise 148 otherwise
149 error('Opperator not supported'); 149 error('Opperator not supported');
150 end 150 end
151 end 151 end
152 152
153 % Closure functions return the opertors applied to the own doamin to close the boundary 153 % Closure functions return the opertors applied to the own doamin to close the boundary
154 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 154 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
155 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 155 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
156 % type is a string specifying the type of boundary condition if there are several. 156 % type is a string specifying the type of boundary condition if there are several.
157 % data is a function returning the data that should be applied at the boundary. 157 % data is a function returning the data that should be applied at the boundary.
165 [closure,penalty] = boundary_condition_general(obj,BM,boundary,L); 165 [closure,penalty] = boundary_condition_general(obj,BM,boundary,L);
166 otherwise 166 otherwise
167 error('No such boundary condition') 167 error('No such boundary condition')
168 end 168 end
169 end 169 end
170 170
171 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 171 function [closure, penalty] = interface(obj, boundary, neighbour_scheme, neighbour_boundary, type)
172 error('An interface function does not exist yet'); 172 error('Not implemented');
173 end 173 end
174 174
175 function N = size(obj) 175 function N = size(obj)
176 N = obj.m; 176 N = obj.m;
177 end 177 end
178 178
179 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y, Z) 179 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y, Z)
180 params = obj.params; 180 params = obj.params;
181 side = max(length(X),length(Y)); 181 side = max(length(X),length(Y));
182 if isa(mat,'function_handle') 182 if isa(mat,'function_handle')
183 [rows,cols] = size(mat(params,0,0,0)); 183 [rows,cols] = size(mat(params,0,0,0));
187 matVec = mat; 187 matVec = mat;
188 [rows,cols] = size(matVec); 188 [rows,cols] = size(matVec);
189 side = max(length(X),length(Y)); 189 side = max(length(X),length(Y));
190 cols = cols/side; 190 cols = cols/side;
191 end 191 end
192 192
193 ret = cell(rows,cols); 193 ret = cell(rows,cols);
194 for ii = 1:rows 194 for ii = 1:rows
195 for jj = 1:cols 195 for jj = 1:cols
196 ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side)); 196 ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side));
197 end 197 end
198 end 198 end
199 ret = cell2mat(ret); 199 ret = cell2mat(ret);
200 end 200 end
201 201
202 function [BM] = boundary_matrices(obj,boundary) 202 function [BM] = boundary_matrices(obj,boundary)
203 params = obj.params; 203 params = obj.params;
204 204
205 switch boundary 205 switch boundary
206 case {'w','W','west'} 206 case {'w','W','west'}
207 BM.e_ = obj.e_w; 207 BM.e_ = obj.e_w;
208 mat = obj.A; 208 mat = obj.A;
209 BM.boundpos = 'l'; 209 BM.boundpos = 'l';
246 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xz,obj.Yz,obj.Z(end)); 246 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.Xz,obj.Yz,obj.Z(end));
247 BM.side = length(obj.Xz); 247 BM.side = length(obj.Xz);
248 end 248 end
249 BM.pos = signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3); 249 BM.pos = signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3);
250 end 250 end
251 251
252 % Characteristic bouyndary consitions 252 % Characteristic bouyndary consitions
253 function [closure, penalty]=boundary_condition_char(obj,BM) 253 function [closure, penalty]=boundary_condition_char(obj,BM)
254 side = BM.side; 254 side = BM.side;
255 pos = BM.pos; 255 pos = BM.pos;
256 neg = BM.neg; 256 neg = BM.neg;
258 V = BM.V; 258 V = BM.V;
259 Vi = BM.Vi; 259 Vi = BM.Vi;
260 Hi = BM.Hi; 260 Hi = BM.Hi;
261 D = BM.D; 261 D = BM.D;
262 e_ = BM.e_; 262 e_ = BM.e_;
263 263
264 switch BM.boundpos 264 switch BM.boundpos
265 case {'l'} 265 case {'l'}
266 tau = sparse(obj.n*side,pos); 266 tau = sparse(obj.n*side,pos);
267 Vi_plus = Vi(1:pos,:); 267 Vi_plus = Vi(1:pos,:);
268 tau(1:pos,:) = -abs(D(1:pos,1:pos)); 268 tau(1:pos,:) = -abs(D(1:pos,1:pos));
274 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); 274 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
275 closure = Hi*e_*V*tau*Vi_minus*e_'; 275 closure = Hi*e_*V*tau*Vi_minus*e_';
276 penalty = -Hi*e_*V*tau*Vi_minus; 276 penalty = -Hi*e_*V*tau*Vi_minus;
277 end 277 end
278 end 278 end
279 279
280 % General boundary condition in the form Lu=g(x) 280 % General boundary condition in the form Lu=g(x)
281 function [closure,penalty] = boundary_condition_general(obj,BM,boundary,L) 281 function [closure,penalty] = boundary_condition_general(obj,BM,boundary,L)
282 side = BM.side; 282 side = BM.side;
283 pos = BM.pos; 283 pos = BM.pos;
284 neg = BM.neg; 284 neg = BM.neg;
285 zeroval=BM.zeroval; 285 zeroval=BM.zeroval;
286 V = BM.V; 286 V = BM.V;
287 Vi = BM.Vi; 287 Vi = BM.Vi;
288 Hi = BM.Hi; 288 Hi = BM.Hi;
289 D = BM.D; 289 D = BM.D;
290 e_ = BM.e_; 290 e_ = BM.e_;
291 291
292 switch boundary 292 switch boundary
293 case {'w','W','west'} 293 case {'w','W','west'}
294 L = obj.evaluateCoefficientMatrix(L,obj.x(1),obj.Yx,obj.Zx); 294 L = obj.evaluateCoefficientMatrix(L,obj.x(1),obj.Yx,obj.Zx);
295 case {'e','E','east'} 295 case {'e','E','east'}
296 L = obj.evaluateCoefficientMatrix(L,obj.x(end),obj.Yx,obj.Zx); 296 L = obj.evaluateCoefficientMatrix(L,obj.x(end),obj.Yx,obj.Zx);
301 case {'b','B','bottom'} 301 case {'b','B','bottom'}
302 L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(1)); 302 L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(1));
303 case {'t','T','top'} 303 case {'t','T','top'}
304 L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(end)); 304 L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(end));
305 end 305 end
306 306
307 switch BM.boundpos 307 switch BM.boundpos
308 case {'l'} 308 case {'l'}
309 tau = sparse(obj.n*side,pos); 309 tau = sparse(obj.n*side,pos);
310 Vi_plus = Vi(1:pos,:); 310 Vi_plus = Vi(1:pos,:);
311 Vi_minus = Vi(pos+zeroval+1:obj.n*side,:); 311 Vi_minus = Vi(pos+zeroval+1:obj.n*side,:);
312 V_plus = V(:,1:pos); 312 V_plus = V(:,1:pos);
313 V_minus = V(:,(pos+zeroval)+1:obj.n*side); 313 V_minus = V(:,(pos+zeroval)+1:obj.n*side);
314 314
315 tau(1:pos,:) = -abs(D(1:pos,1:pos)); 315 tau(1:pos,:) = -abs(D(1:pos,1:pos));
316 R = -inv(L*V_plus)*(L*V_minus); 316 R = -inv(L*V_plus)*(L*V_minus);
317 closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; 317 closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
318 penalty = -Hi*e_*V*tau*inv(L*V_plus)*L; 318 penalty = -Hi*e_*V*tau*inv(L*V_plus)*L;
319 case {'r'} 319 case {'r'}
320 tau = sparse(obj.n*side,neg); 320 tau = sparse(obj.n*side,neg);
321 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); 321 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
322 Vi_plus = Vi(1:pos,:); 322 Vi_plus = Vi(1:pos,:);
323 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:); 323 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
324 324
325 V_plus = V(:,1:pos); 325 V_plus = V(:,1:pos);
326 V_minus = V(:,(pos+zeroval)+1:obj.n*side); 326 V_minus = V(:,(pos+zeroval)+1:obj.n*side);
327 R = -inv(L*V_minus)*(L*V_plus); 327 R = -inv(L*V_minus)*(L*V_plus);
328 closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_'; 328 closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
329 penalty = -Hi*e_*V*tau*inv(L*V_minus)*L; 329 penalty = -Hi*e_*V*tau*inv(L*V_minus)*L;
330 end 330 end
331 end 331 end
332 332
333 % Function that diagonalizes a symbolic matrix A as A=V*D*Vi 333 % Function that diagonalizes a symbolic matrix A as A=V*D*Vi
334 % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign 334 % D is a diagonal matrix with the eigenvalues on A on the diagonal sorted by their sign
335 % [d+ ] 335 % [d+ ]
336 % D = [ d0 ] 336 % D = [ d0 ]
337 % [ d-] 337 % [ d-]
342 [V, D] = eig(mat(params,xs,ys,zs)); 342 [V, D] = eig(mat(params,xs,ys,zs));
343 Vi=inv(V); 343 Vi=inv(V);
344 xs = x; 344 xs = x;
345 ys = y; 345 ys = y;
346 zs = z; 346 zs = z;
347 347
348 348
349 side = max(length(x),length(y)); 349 side = max(length(x),length(y));
350 Dret = zeros(obj.n,side*obj.n); 350 Dret = zeros(obj.n,side*obj.n);
351 Vret = zeros(obj.n,side*obj.n); 351 Vret = zeros(obj.n,side*obj.n);
352 Viret= zeros(obj.n,side*obj.n); 352 Viret= zeros(obj.n,side*obj.n);
353 353
354 for ii=1:obj.n 354 for ii=1:obj.n
355 for jj=1:obj.n 355 for jj=1:obj.n
356 Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii)); 356 Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii));
357 Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii)); 357 Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii));
358 Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii)); 358 Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii));
359 end 359 end
360 end 360 end
361 361
362 D = sparse(Dret); 362 D = sparse(Dret);
363 V = sparse(Vret); 363 V = sparse(Vret);
364 Vi = sparse(Viret); 364 Vi = sparse(Viret);
365 V = obj.evaluateCoefficientMatrix(V,x,y,z); 365 V = obj.evaluateCoefficientMatrix(V,x,y,z);
366 Vi= obj.evaluateCoefficientMatrix(Vi,x,y,z); 366 Vi= obj.evaluateCoefficientMatrix(Vi,x,y,z);
367 D = obj.evaluateCoefficientMatrix(D,x,y,z); 367 D = obj.evaluateCoefficientMatrix(D,x,y,z);
368 DD = diag(D); 368 DD = diag(D);
369 369
370 poseig = (DD>0); 370 poseig = (DD>0);
371 zeroeig = (DD==0); 371 zeroeig = (DD==0);
372 negeig = (DD<0); 372 negeig = (DD<0);
373 373
374 D = diag([DD(poseig); DD(zeroeig); DD(negeig)]); 374 D = diag([DD(poseig); DD(zeroeig); DD(negeig)]);
375 V = [V(:,poseig) V(:,zeroeig) V(:,negeig)]; 375 V = [V(:,poseig) V(:,zeroeig) V(:,negeig)];
376 Vi= [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)]; 376 Vi= [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
377 signVec = [sum(poseig),sum(zeroeig),sum(negeig)]; 377 signVec = [sum(poseig),sum(zeroeig),sum(negeig)];
378 end 378 end