LGL_EIG Computes eigenvalues of first/second order spectral derivative matrices: collocation/G-NI, LGL nodes, 1D monodomain [d,A]=lgl_eig(nx,nu,pbl); Input: nx = polynomial degree nu = viscosity (constant > 0) pbl = integer parameter to set the problem: 0 : compute eigenvalues of A with Au := du/dx in (-1,1), u(1)=0 strong form (collocation) = 1 : compute eigenvalues of A with Au := du/dx in (-1,1), u(1)=0 weak form with integration by part = 2 : compute generalized eigenvalues of A with Au := du/dx in (-1,1), u(1)=0 weak form with integration by part = 30 : compute eigenvalues of A with Au := -d^2u/dx^2 in (-1,1), u(-1)=u(1)=0 strong form = 31 : compute eigenvalues of A with Au := -d^2u/dx^2 in (-1,1), u(-1)=u(1)=0 weak form (G-NI) = 32 : compute generalized eigenvalues of A with Au := -d^2u/dx^2 in (-1,1), u(-1)=u(1)=0 weak form (G-NI) = 41 : compute eigenvalues of A with Au := -d^2u/dx^2 in (-1,1), u'(-1)=u'(1)=0 weak form (G-NI) = 42 : compute generalized eigenvalues of A with Au := -d^2u/dx^2 in (-1,1), u'(-1)=u'(1)=0 weak form (G-NI) = 51 : compute eigenvalues of A with Au := -nu d^2u/dx^2+du/dx in (-1,1), u(-1)=u(1)=0 weak form (G-NI) = 52 : compute generalized eigenvalues of A with Au := -nu d^2u/dx^2+du/dx in (-1,1), u(-1)=u(1)=0 weak form (G-NI) = 61 : compute eigenvalues of A with Au := -nu d^2u/dx^2+du/dx in (-1,1), u(-1)=0, u'(1)=0 weak form (G-NI) = 62 : compute generalized eigenvalues of A with Au := -nu d^2u/dx^2+du/dx in (-1,1), u(-1)=0, u'(1)=0 weak form (G-NI) Output: d = column array with eigenvalues of A A = matrix required, of size (npdx,npdx) Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, "Spectral Methods. Fundamentals in Single Domains" Springer Verlag, Berlin Heidelberg New York, 2006.
0001 function [d,A]=lgl_eig(nx,nu,pbl); 0002 % LGL_EIG Computes eigenvalues of first/second order spectral derivative matrices: collocation/G-NI, LGL nodes, 1D monodomain 0003 % 0004 % [d,A]=lgl_eig(nx,nu,pbl); 0005 % 0006 % Input: nx = polynomial degree 0007 % nu = viscosity (constant > 0) 0008 % pbl = integer parameter to set the problem: 0009 % 0 : compute eigenvalues of A with 0010 % Au := du/dx in (-1,1), u(1)=0 0011 % strong form (collocation) 0012 % 0013 % = 1 : compute eigenvalues of A with 0014 % Au := du/dx in (-1,1), u(1)=0 0015 % weak form with integration by part 0016 % 0017 % = 2 : compute generalized eigenvalues of A with 0018 % Au := du/dx in (-1,1), u(1)=0 0019 % weak form with integration by part 0020 % 0021 % = 30 : compute eigenvalues of A with 0022 % Au := -d^2u/dx^2 in (-1,1), u(-1)=u(1)=0 0023 % strong form 0024 % 0025 % = 31 : compute eigenvalues of A with 0026 % Au := -d^2u/dx^2 in (-1,1), u(-1)=u(1)=0 0027 % weak form (G-NI) 0028 % 0029 % = 32 : compute generalized eigenvalues of A with 0030 % Au := -d^2u/dx^2 in (-1,1), u(-1)=u(1)=0 0031 % weak form (G-NI) 0032 % 0033 % = 41 : compute eigenvalues of A with 0034 % Au := -d^2u/dx^2 in (-1,1), u'(-1)=u'(1)=0 0035 % weak form (G-NI) 0036 % 0037 % = 42 : compute generalized eigenvalues of A with 0038 % Au := -d^2u/dx^2 in (-1,1), u'(-1)=u'(1)=0 0039 % weak form (G-NI) 0040 % 0041 % = 51 : compute eigenvalues of A with 0042 % Au := -nu d^2u/dx^2+du/dx in (-1,1), u(-1)=u(1)=0 0043 % weak form (G-NI) 0044 % 0045 % = 52 : compute generalized eigenvalues of A with 0046 % Au := -nu d^2u/dx^2+du/dx in (-1,1), u(-1)=u(1)=0 0047 % weak form (G-NI) 0048 % 0049 % = 61 : compute eigenvalues of A with 0050 % Au := -nu d^2u/dx^2+du/dx in (-1,1), u(-1)=0, u'(1)=0 0051 % weak form (G-NI) 0052 % 0053 % = 62 : compute generalized eigenvalues of A with 0054 % Au := -nu d^2u/dx^2+du/dx in (-1,1), u(-1)=0, u'(1)=0 0055 % weak form (G-NI) 0056 % 0057 % Output: d = column array with eigenvalues of A 0058 % A = matrix required, of size (npdx,npdx) 0059 % 0060 % Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0061 % "Spectral Methods. Fundamentals in Single Domains" 0062 % Springer Verlag, Berlin Heidelberg New York, 2006. 0063 0064 % Written by Paola Gervasio 0065 % $Date: 2007/04/01$ 0066 0067 % 0068 npdx=nx+1; 0069 [x,wx]=xwlgl(npdx); dx=derlgl(x,npdx); 0070 0071 % A 0072 0073 % 0 : compute eigenvalues of A with 0074 % Au := du/dx in (-1,1), u(1)=0 0075 % strong form (collocation) 0076 if pbl == 0 0077 A=dx(1:nx,1:nx); 0078 0079 elseif pbl == 1 0080 % = 1 : compute eigenvalues of A with 0081 % Au := du/dx in (-1,1), u(1)=0 0082 % weak form with integration by part 0083 A=-dx'*diag(wx); 0084 A(1,1)=A(1,1)-1; 0085 0086 elseif pbl == 2 0087 % = 2 : compute generalized eigenvalues of A with 0088 % Au := du/dx in (-1,1), u(1)=0 0089 % weak form with integration by part 0090 A=-dx'*diag(wx); 0091 A(1,1)=A(1,1)-1; 0092 A=diag(1./wx)*A; 0093 0094 elseif pbl == 30 0095 % = 30 : compute eigenvalues of A with 0096 % Au := -d^2u/dx^2 in (-1,1), u(-1)=u(1)=0 0097 % strong form 0098 A=-dx*dx; 0099 A=A(2:nx,2:nx); 0100 0101 elseif pbl == 31 0102 % = 31 : compute eigenvalues of A with 0103 % Au := -d^2u/dx^2 in (-1,1), u(-1)=u(1)=0 0104 % weak form (G-NI) 0105 A=dx'*diag(wx)*dx; 0106 A=A(2:nx,2:nx); 0107 0108 elseif pbl == 32 0109 % = 32 : compute generalized eigenvalues of A with 0110 % Au := -d^2u/dx^2 in (-1,1), u(-1)=u(1)=0 0111 % weak form (G-NI) 0112 A=dx'*diag(wx)*dx; 0113 A=diag(1./wx(2:nx))*A(2:nx,2:nx); 0114 d=eig(A); 0115 0116 elseif pbl == 41 0117 % = 41 : compute eigenvalues of A with 0118 % Au := -d^2u/dx^2 in (-1,1), u'(-1)=u'(1)=0 0119 % weak form (G-NI) 0120 A=dx'*diag(wx)*dx; 0121 A=A(1:npdx,1:npdx); 0122 0123 elseif pbl == 42 0124 % = 42 : compute generalized eigenvalues of A with 0125 % Au := d^2u/dx^2 in (-1,1), u'(-1)=u'(1)=0 0126 % weak form (G-NI) 0127 % 0128 A=dx'*diag(wx)*dx; 0129 A=diag(1./wx)*A; 0130 0131 elseif pbl == 51 0132 % = 51 : compute eigenvalues of A with 0133 % Au := -nu d^2u/dx^2+du/dx in (-1,1), u(-1)=u(1)=0 0134 % weak form (G-NI) 0135 A=dx'*diag(wx)*dx; 0136 A=A(2:nx,2:nx); 0137 0138 elseif pbl == 52 0139 % = 52 : compute generalized eigenvalues of A with 0140 % Au := -nu d^2u/dx^2+du/dx in (-1,1), u(-1)=u(1)=0 0141 % weak form (G-NI) 0142 A=dx'*diag(wx)*dx;A=A*nu; 0143 A=A-dx'*diag(wx); 0144 A=diag(1./wx)*A; 0145 A=A(2:nx,2:nx); 0146 0147 elseif pbl == 61 0148 % = 61 : compute eigenvalues of A with 0149 % Au := -nu d^2u/dx^2+du/dx in (-1,1), u(-1)=0, u'(1)=0 0150 % weak form (G-NI) 0151 A=dx'*diag(wx)*dx;A=A*nu; 0152 A=A-dx'*diag(wx); 0153 A(npdx,npdx)=A(npdx,npdx)+1; 0154 A=A(2:npdx,2:npdx); 0155 0156 elseif pbl == 62 0157 % = 62 : compute generalized eigenvalues of A with 0158 % Au := -nu d^2u/dx^2+du/dx in (-1,1), u(-1)=0, u'(1)=0 0159 % weak form (G-NI) 0160 A=dx'*diag(wx)*dx;A=A*nu; 0161 A=A-dx'*diag(wx); 0162 A(npdx,npdx)=A(npdx,npdx)+1; 0163 A=diag(1./wx)*A; 0164 A=A(2:npdx,2:npdx); 0165 end 0166 d=eig(A);