Spectral Integral Suite in C++
testRiblets.cpp
Go to the documentation of this file.
1 
2 #include <fstream> // To output data to files
3 #include <iostream>
4 #include <stdio.h>
5 #include <sis.hpp>
6 
7 #include <eigen3/unsupported/Eigen/MatrixFunctions>
8 
9 using namespace std;
10 #ifndef CD_T
11 typedef complex<double> Cd_t;
12 #endif
13 
14 typedef valarray<Cd_t> Vcd_t;
15 //#define DEBUG
16 
17 #include <lyap.h>
18 
19 int main()
20 {
21  //Tau formulation, discriptor system :
22  int N = 91;
23  Cd_t ii(0,1); // imaginary unit.
24  Eigen::VectorXd kzval = Eigen::VectorXd::LinSpaced(10, 1.0, 4.0);
25  // for (int i = 0; i < 10; i++){
26  double kx = 1;
27  //double kz = kzval(i);
28  double kz = 1.0;
29  double k2 = kx * kx + kz * kz;
30  double k4 = k2 * k2;
31  double R = 2000;
32  Eigen::MatrixXcd Dmat = Eigen::MatrixXcd::Zero(N + 1, N + 1);
33  Eigen::MatrixXcd Dpmat = Eigen::MatrixXcd::Zero(N - 1, N - 1);
34  for (int k = 0; k < N + 1; k++) {
35  for (int p = k + 1; p < N + 1; p++) {
36  if ((p + k) % 2) {
37  Dmat(k, p) = p;
38  }
39  }
40  }
41  cout << "line: " << __LINE__ << endl << flush;
42  for (int k = 0; k < N - 1; k++) {
43  for (int p = k + 1; p < N - 1; p++) {
44  if ((p + k) % 2) {
45  Dpmat(k, p) = p;
46  }
47  }
48  }
49  cout << "line: " << __LINE__ << endl
50  << flush;
51  Dmat = Dmat * 2.0;
52  Dpmat = Dpmat * 2.0;
53  Eigen::MatrixXcd If = Eigen::MatrixXcd::Identity(N + 1, N + 1);
54  Eigen::MatrixXcd Ip = Eigen::MatrixXcd::Identity(N - 1, N - 1);
55  Eigen::MatrixXcd D1f = Dmat * If;
56  Eigen::MatrixXcd D2f = D1f * D1f;
57  Eigen::MatrixXcd D1p = Dpmat * Ip;
58  Eigen::MatrixXcd D2p = Dpmat * Dpmat;
59  cout << "line: " << __LINE__ << endl
60  << flush;
61  // Make discretization matrices : % First verify eigenvalues :
62  Eigen::MatrixXcd bc1(1, N + 1);
63  bc1.setConstant(1.0);
64  Eigen::MatrixXcd bc2(1, N + 1);
65  for (int i = 0; i < N + 1; i++){
66  bc2(0, i) = pow(-1.0, i);
67  }
68  //bc2 = (-1).^ (0: N);
69  bc1(0,0) = 0.5 * bc1(0,0);
70  bc2(0,0) = 0.5 * bc2(0,0);
71  cout << "line: " << __LINE__ << endl
72  << flush;
73  Eigen::MatrixXcd I = If.block(0, 0, N - 1, N + 1);
74  //I = If(1: N - 1,:);
75  Eigen::MatrixXcd D1 = D1f.block(0, 0, N - 1, N + 1);
76  //D1 = D1f(1: N - 1,:);
77  Eigen::MatrixXcd D2 = D2f.block(0, 0, N - 1, N + 1);
78  //D2 = D2f(1: N - 1,:);
79  Eigen::MatrixXcd E(4 * (N - 1) + 6, 4 * (N - 1) + 6), F(4 * (N - 1) + 6, 4 * (N - 1) + 6), B(4 * (N - 1) + 6, 3 * (N - 1)), C(3 * (N + 1), 3 * (N + 1) + N - 1);
80  cout << "line: " << __LINE__ << endl
81  << flush;
82  E << I, 0.0 * I, I * 0.0, 0.0 * Ip, //
83  0.0 * I, I, I * 0.0, 0.0 * Ip, //
84  0.0 * I, 0.0 * I, I , 0.0 * Ip, //
85  0.0 * I, 0.0 * I, I * 0.0, 0.0 * Ip, //
86  Eigen::MatrixXcd::Zero(6, 3 * (N + 1) + N - 1);
87  // E = [ I, 0 * I, 0 * I, 0 * Ip;
88  // 0 * I, I, 0 * I, 0 * Ip;
89  // 0 * I, 0 * I, I, 0 * Ip;
90  // 0 * I, 0 * I, 0 * I, 0 * Ip;
91  // zeros(6, 3 * (N + 1) + N - 1) ];
92  cout << "line: " << __LINE__ << endl
93  << flush;
94  sis::N = N-2;
95  sis::sis_setup();
96  cout << "line: " << __LINE__ << endl
97  << flush;
98  Eigen::MatrixXcd U = sis::Chebfun<Cd_t>(Vcd_t(1.0 - pow(sis::yc, 2))).MultMat();
99  Eigen::MatrixXcd Uy = -sis::Chebfun<Cd_t>(Vcd_t(2.0 * sis::yc)).MultMat();
100  cout << "line: " << __LINE__ << endl
101  << flush;
102  // U = MultMat(1 - y2.^ 2);
103  // Uy = -MultMat(2 * y2);
104  Eigen::MatrixXcd Delta = D2 - k2 * I;
105  Eigen::MatrixXcd Diag = (1 / R) * Delta - ii * kx * U * I;
106  cout << "line: " << __LINE__ << endl
107  << flush;
108  F << Diag, -Uy * I, I * 0, -ii * kx * Ip,//
109  I * 0, Diag, I * 0, -D1p,//
110  I * 0, I * 0, Diag, -ii * kz * Ip,//
111  ii * kx * I, D1, ii * kz * I, 0 * Ip,//
112  bc1, bc1 * 0, bc1 * 0, Eigen::MatrixXcd::Zero(1, N - 1),//
113  bc2, bc2 * 0, bc2 * 0, Eigen::MatrixXcd::Zero(1, N - 1),//
114  0 * bc1, bc1, bc1 * 0, Eigen::MatrixXcd::Zero(1, N - 1),//
115  0 * bc2, bc2, bc2 * 0, Eigen::MatrixXcd::Zero(1, N - 1),//
116  0 * bc1, bc1 * 0, bc1, Eigen::MatrixXcd::Zero(1, N - 1),//
117  0 * bc2, bc2 * 0, bc2, Eigen::MatrixXcd::Zero(1, N - 1);
118  cout << "line: " << __LINE__ << endl
119  << flush;
120  B << Ip, Ip*0, Ip*0,//
121  Ip*0,Ip, Ip*0,//
122  Ip*0, Ip*0, Ip,//
123  Eigen::MatrixXcd::Zero(N-1+6, 3*(N-1));
124  cout << "line: " << __LINE__ << endl
125  << flush;
126  // B = [blkdiag(eye(N - 1), eye(N - 1), eye(N - 1)); zeros((N - 1) + 6, 3 * (N - 1))];
127 
128  C << If, 0 * If, 0 * If, Eigen::MatrixXcd::Zero(N + 1, N - 1), //
129  0 * If, If, 0 * If, Eigen::MatrixXcd::Zero(N + 1, N - 1), //
130  0 * If, 0 * If, If, Eigen::MatrixXcd::Zero(N + 1, N - 1);
131  //C = [ If, 0 * If, 0 * If, zeros(N + 1, N - 1);
132  // 0 * If, If, 0 * If, zeros(N + 1, N - 1);
133  // 0 * If, 0 * If, If, zeros(N + 1, N - 1) ];
134  cout << "line: " << __LINE__ << endl
135  << flush;
136  // Verify eigensystem:
137  OrdQz qz(F,E);
138  cout << "line: " << __LINE__ << endl
139  << flush;
140  Eigen::MatrixXcd evals = qz.alpha.array()/qz.beta.array();
141  cout << "line: " << __LINE__ << endl
142  << flush;
143  cout << "Showing the eigenvalues: " << evals << endl;
144  int bre;
145  cin >> bre;
146  cout << "info: " << qz.info << endl;
147  Eigen::ColPivHouseholderQR<Eigen::Matrix<std::complex<double>, Eigen::Dynamic, Eigen::Dynamic> > qr;
148  qr.compute(qz.T.block(0, 0, qz.M, qz.M));
149  cout << "is invertible: " << qr.isInvertible() << endl;
150  //for (int i = 0; i < qz.alpha.size(); i++){
151  // cout << criteria_(qz.alpha(i,0),qz.beta(i,0)) << endl;
152  // cout << "qz.beta: " << qz.beta(i,0) << endl;
153  //}
154 /*
155  Eigen::MatrixXcd Iw = intWts(N);
156  Eigen::MatrixXcd Iwp = intWts(N - 2);
157 
158  S1p = blkdiag(Iwp, Iwp, Iwp);
159  sqrtmS1p = sqrtm(S1p);
160  invsqrtmS1p = sqrtmS1p\eye(size(sqrtmS1p));
161 
162  S1 = blkdiag(Iw, Iw, Iw);
163  sqrtmS1 = sqrtm(S1);
164  invsqrtmS1 = sqrtmS1\eye(size(sqrtmS1));
165 
166  S2 = blkdiag(Iw, Iw, Iw, Iwp);
167  sqrtmS2 = sqrtm(S2);
168  invsqrtmS2 = sqrtmS2\eye(size(sqrtmS2));
169  invIw = eye(N + 1)\Iw;
170 
171  % In new coordinates : Es = sqrtmS2 *E *invsqrtmS2;
172  Fs = sqrtmS2 * F * invsqrtmS2;
173  Bs = sqrtmS2 * B * invsqrtmS1p;
174  Cs = sqrtmS1 * C * invsqrtmS2;
175  Es = complex(Es, zeros(size(Es)));
176 
177  [ AAS, BBS, Qs, Zs, alp, bet, info ] = ordqzLapack(Fs, Es);
178 
179  % [ AA, BB, Q, Z, V, W ] = qz(Fs, Es);
180  % [ AAS, BBS, QS, ZS ] = ordqz(AA, BB, Q, Z, 'lhp');
181  evals = alp./ bet;
182  [ r, c ] = size(AAS);
183  isinfE = isinf(evals);
184  sumisinfE = sum(isinfE);
185  p = r - sumisinfE;
186  W = QS '; T = ZS';
187  Ff = AAS(1
188  : p, 1
189  : p);
190  Fu = AAS(1
191  : p, p + 1
192  : end);
193  Finf = AAS(p + 1
194  : end, p + 1
195  : end);
196 
197  Ef = BBS(1
198  : p, 1
199  : p);
200  Eu = BBS(1
201  : p, p + 1
202  : end);
203  Einf = BBS(p + 1
204  : end, p + 1
205  : end);
206  [ Y, Z ] = solveCoupledSylvester(Ef, Einf, -Eu, Ff, Finf, -Fu);
207  % check norms : norm(Ef * Y - Z * Einf + Eu)
208  norm(Ff * Y - Z * Finf + Fu)
209  [r, c] = size(Z);
210  [ ry, cy ] = size(Y);
211  [ r1, c1 ] = size(W);
212  [ re, ce ] = size(Es);
213  [ rt, ct ] = size(T);
214  V = W * [ eye(r, re - c), Z; zeros(r1 - r, re - c), eye(r1 - r, c) ];
215  U = [ eye(ry, rt - cy), -Y; zeros(ce - ry, r1 - cy), eye(ce - ry, c) ] * T;
216  % Check norms : norm(Es - V * blkdiag(Ef, Einf) * U)
217  norm(Fs - V * blkdiag(Ff, Finf) * U)
218 
219  [rv, cv] = size(V);
220  [ ref, cef ] = size(Ef);
221  [ rfinf, cfinf ] = size(Finf);
222 
223  % Solve generalized Lyapunov equations : invEfFf = Ef\Ff;
224 Ftilde = W'*(Bs*Bs')*W;
225 Ftilde1 = Ftilde(1
226  : ref, 1
227  : cef);
228 Ftilde2 = Ftilde(1
229  : ref, cef + 1
230  : end);
231 Ftilde3 = Ftilde(ref + 1
232  : end, 1
233  : cef);
234 Ftilde4 = Ftilde(ref + 1
235  : end, cef + 1
236  : end);
237 Gc1 = lyap(invEfFf, Ftilde1 - Ftilde2 * Z ' - Z*Ftilde3 + Z*Ftilde4*Z');
238 Gc1 = Ef\(Gc1/(Ef'));
239 Gc = T'*[Gc1, zeros(ref,ct-cef);
240  zeros(ct-cef,cef),zeros(ct-cef,ct-cef)]*T;
241 invFinfEinf = Finf\Einf;
242 Gnc4 = dlyap(invFinfEinf,-Ftilde4);
243 Gnc4 = Finf\(Gnc4/(Finf'));
244 Gnc = T'*[Y*Gnc4*Y', Y*Gnc4;
245  Gnc4*Y', Gnc4]*T;
246 
247 Gc = sqrtmS2*Gc*invsqrtmS2;
248 Gnc = sqrtmS2*Gnc*invsqrtmS2;
249 
250  normVec(i) = trace(Cs*(Gc+Gnc)*Cs');
251 end
252 
253 
254 
255 % Iw4 = blkdiag(Iw,Iw,Iw,Iw);
256 % Iw3 = blkdiag(Iw,Iw,Iw);
257 % invIw4 = eye((N+1)*4)\Iw4;
258 % invIw3 = eye((N+1)*3)\Iw3;
259 % nullbc1 = blkdiag(nullbcs,nullbcs,nullbcs,eye(90,90));
260 % % Now come to calculating the H2 norm:
261 % C = [If, 0*If, 0*If, zeros(92,90);
262 % 0*If, If, 0*If, zeros(92,90);
263 % 0*If, 0*If, If, zeros(92,90)]*nullbc1*nullcont;
264 % Q = C'*Iw3*C;
265 
266 */
267 
268 
269  return 0;
270 }
sis::sis_setup
void sis_setup()
Definition: sis.hpp:954
OrdQz::beta
Eigen::MatrixXcd beta
Definition: lyap.h:36
sis::Chebfun
This is a chebfun analogue. Chebfun will represent both values in physical space or an array of Cheby...
Definition: sis.hpp:1150
OrdQz
Definition: lyap.h:34
std
Definition: sis.hpp:260
ii
complex< double > ii(0.0, 1.0)
OrdQz::info
int info
Definition: lyap.h:37
OrdQz::T
Eigen::MatrixXcd T
Definition: lyap.h:36
sis::yc
std::valarray< std::complex< SIS_TYPE > > yc(N+1)
Cd_t
complex< double > Cd_t
Definition: testRiblets.cpp:11
std::pow
valarray< complex< T > > pow(valarray< complex< T > > base, T power)
Definition: sis.hpp:453
sis::N
int N
Specifies number of Chebyshev polynomials, default N = 31.
Definition: sis.hpp:472
main
int main()
Definition: testRiblets.cpp:19
OrdQz::M
int M
Definition: lyap.h:38
Vcd_t
valarray< Cd_t > Vcd_t
Definition: testRiblets.cpp:14
OrdQz::alpha
Eigen::MatrixXcd alpha
Definition: lyap.h:36
Cd_t
complex< double > Cd_t
Definition: Ex_05.cpp:17