ELLPRECOFEM_1D FEM-preconditioned SEM-appx of 1D second order b.v.p.
Numerical solution (or eigenvalues compuation)
of the 1D boundary value problem
-nu u''+ beta*u' + gam u =f xa < x < xb
Dirichlet or Neumann boundary conditions
by Galerkin Numerical Integration with LGL quadrature formulas
preconditioner by FEM matrices
Useful to reproduce results on iterative condition numbers (1D case)
published in Ch4. of CHQZ2.
[xy,un,A,M,AFE,MFE,MFEd,d,kappa,param]=ellprecofem_1d(xa,xb,nu,...
beta, gam,ff,cb,ub,ne,nx,param);
Input: xa, xb = extrema of computational domain Omega=(xa,xb)
nu = viscosity (constant>0)
beta = coefficient of first order term (constant)
gam = coefficient of zeroth order term (constant>=0)
ff = r.h.s. solution (ff=@(x)[ff(x)], with .*, .^, ./)
cb = string containing boundary conditions, for example
cb='dn' imposes Dirichlet in xa and Neumann in xb,
possible b.c.: Dirichlet, Neumann
ub = 2 components array with Dirichlet, Neumann (pure derivative)
data in xa and in xb
ne = number of elements (equally spaced)
nx = polynomial degree in each element (the same in each element)
to be set only if p=4, otherwise, nx=p;
param = parameters array :
param(1)= choice of the preconditioner:
0 : P=I
A=K_GNI
1 : P=K_FE (4.4.45), pag. 221, CHQZ2
A=K_GNI
2 : P=M_FE^{-1}K_FE (4.4.46), pag. 221, CHQZ2
A=M_GNI^{-1}K_GNI
3 : P=M_FE,d^{-1}K_FE (4.4.47), pag. 221, CHQZ2
A=M_GNI^{-1}K_GNI
4 : P=M_FE^{-1/2}K_FE M_FE^{-1/2} (4.4.48), pag. 221, CHQZ2
A=M_GNI^{-1/2}K_GNI M_GNI^{-1/2}
5 : P=KM_FE<d^{-1/2}K_FE M_FE,d^{-1/2}_FE (4.4.49), pag. 221, CHQZ2
A=M_GNI^{-1/2}K_GNI M_GNI^{-1/2}
param(2)= choice of the solver:
1 : Conjugate Gradient, only for symmetric preconditioners:
usable only if param(1)=1, 4, 5
2 : Bi-CGStab (VDV), usable for any choice of param(1)
11 : eigenvalues computation, to have iterative condition
number
kapit(P^{-1}A)=|lambda_max(P^{-1}A)|/|lambda_min(P^{-1}A)|
12 : singular values computation, to have 2-norm condition
number
k_2(P^{-1}A)=sigma_max(P^{-1}A)/sigma_min(P^{-1}A)
param(3) = tolerance for the solver:
stopping test is: ||r^k||_2/||r^0||_2 < tol
(not used if param(2)=11 | param(2)=12)
param(4) = maximum iteration number for the solver
(not used if param(2)=11 | param(2)=12)
Output: xy = 1D mesh
un = numerical solution
A = Stiffness matrix K_GNI (formula (4.4.43), pag. 219 CHQZ2)
M = Mass matrix K_GNI (formula (4.4.43), pag. 219 CHQZ2)
AFE = FEM Stiffness matrix K_FE (Tab. 4.6, pag. 221, CHQZ2)
MFE = FEM Mass matrix (Tab. 4.6, pag. 221, CHQZ2)
MFE = FEM Mass matrix with numerical integration
(Tab. 4.6, pag. 221, CHQZ2)
d = eigenvalues or singular values of the preconditioned matrix
kappa = iterative or 2-norm condition number (if param(2)=11, 12)
param(11) = iterations required by the solver
param(12) = residual at the last iterations
References: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang,
"Spectral Methods. Fundamentals in Single Domains"
Springer Verlag, Berlin Heidelberg New York, 2006.
VDV = van der Vorst, Henk A.,
"Iterative Krylov methods for large linear systems"
Cambridge Monographs on Applied and Computational Mathematics, 13.
Cambridge University Press, Cambridge, 2003.0001 function [xy,un,A,M,AFE,MFE,MFEd,d,kappa,param]=ellprecofem_1d(xa,xb,nu,... 0002 beta, gam,ff,cb,ub,ne,nx,param); 0003 % ELLPRECOFEM_1D FEM-preconditioned SEM-appx of 1D second order b.v.p. 0004 % 0005 % Numerical solution (or eigenvalues compuation) 0006 % of the 1D boundary value problem 0007 % 0008 % -nu u''+ beta*u' + gam u =f xa < x < xb 0009 % Dirichlet or Neumann boundary conditions 0010 % 0011 % by Galerkin Numerical Integration with LGL quadrature formulas 0012 % preconditioner by FEM matrices 0013 % 0014 % Useful to reproduce results on iterative condition numbers (1D case) 0015 % published in Ch4. of CHQZ2. 0016 % 0017 % 0018 % [xy,un,A,M,AFE,MFE,MFEd,d,kappa,param]=ellprecofem_1d(xa,xb,nu,... 0019 % beta, gam,ff,cb,ub,ne,nx,param); 0020 % 0021 % Input: xa, xb = extrema of computational domain Omega=(xa,xb) 0022 % nu = viscosity (constant>0) 0023 % beta = coefficient of first order term (constant) 0024 % gam = coefficient of zeroth order term (constant>=0) 0025 % ff = r.h.s. solution (ff=@(x)[ff(x)], with .*, .^, ./) 0026 % cb = string containing boundary conditions, for example 0027 % cb='dn' imposes Dirichlet in xa and Neumann in xb, 0028 % possible b.c.: Dirichlet, Neumann 0029 % ub = 2 components array with Dirichlet, Neumann (pure derivative) 0030 % data in xa and in xb 0031 % ne = number of elements (equally spaced) 0032 % nx = polynomial degree in each element (the same in each element) 0033 % to be set only if p=4, otherwise, nx=p; 0034 % 0035 % param = parameters array : 0036 % param(1)= choice of the preconditioner: 0037 % 0 : P=I 0038 % A=K_GNI 0039 % 1 : P=K_FE (4.4.45), pag. 221, CHQZ2 0040 % A=K_GNI 0041 % 2 : P=M_FE^{-1}K_FE (4.4.46), pag. 221, CHQZ2 0042 % A=M_GNI^{-1}K_GNI 0043 % 3 : P=M_FE,d^{-1}K_FE (4.4.47), pag. 221, CHQZ2 0044 % A=M_GNI^{-1}K_GNI 0045 % 4 : P=M_FE^{-1/2}K_FE M_FE^{-1/2} (4.4.48), pag. 221, CHQZ2 0046 % A=M_GNI^{-1/2}K_GNI M_GNI^{-1/2} 0047 % 5 : P=KM_FE<d^{-1/2}K_FE M_FE,d^{-1/2}_FE (4.4.49), pag. 221, CHQZ2 0048 % A=M_GNI^{-1/2}K_GNI M_GNI^{-1/2} 0049 % param(2)= choice of the solver: 0050 % 1 : Conjugate Gradient, only for symmetric preconditioners: 0051 % usable only if param(1)=1, 4, 5 0052 % 2 : Bi-CGStab (VDV), usable for any choice of param(1) 0053 % 11 : eigenvalues computation, to have iterative condition 0054 % number 0055 % kapit(P^{-1}A)=|lambda_max(P^{-1}A)|/|lambda_min(P^{-1}A)| 0056 % 12 : singular values computation, to have 2-norm condition 0057 % number 0058 % k_2(P^{-1}A)=sigma_max(P^{-1}A)/sigma_min(P^{-1}A) 0059 % param(3) = tolerance for the solver: 0060 % stopping test is: ||r^k||_2/||r^0||_2 < tol 0061 % (not used if param(2)=11 | param(2)=12) 0062 % param(4) = maximum iteration number for the solver 0063 % (not used if param(2)=11 | param(2)=12) 0064 % 0065 % 0066 % Output: xy = 1D mesh 0067 % un = numerical solution 0068 % A = Stiffness matrix K_GNI (formula (4.4.43), pag. 219 CHQZ2) 0069 % M = Mass matrix K_GNI (formula (4.4.43), pag. 219 CHQZ2) 0070 % AFE = FEM Stiffness matrix K_FE (Tab. 4.6, pag. 221, CHQZ2) 0071 % MFE = FEM Mass matrix (Tab. 4.6, pag. 221, CHQZ2) 0072 % MFE = FEM Mass matrix with numerical integration 0073 % (Tab. 4.6, pag. 221, CHQZ2) 0074 % d = eigenvalues or singular values of the preconditioned matrix 0075 % kappa = iterative or 2-norm condition number (if param(2)=11, 12) 0076 % param(11) = iterations required by the solver 0077 % param(12) = residual at the last iterations 0078 % 0079 % 0080 % References: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0081 % "Spectral Methods. Fundamentals in Single Domains" 0082 % Springer Verlag, Berlin Heidelberg New York, 2006. 0083 % VDV = van der Vorst, Henk A., 0084 % "Iterative Krylov methods for large linear systems" 0085 % Cambridge Monographs on Applied and Computational Mathematics, 13. 0086 % Cambridge University Press, Cambridge, 2003. 0087 0088 % Written by Paola Gervasio 0089 % $Date: 2007/04/01$ 0090 0091 % check if the choice for the preconditioner and the solver are compatible 0092 xy=[];un=[];A=[];M=[];AFE=[];MFE=[];MFEd=[];d=[];kappa=1; 0093 0094 if param(2)==1 0095 if(param(1)==2 | param(1)==3) 0096 fprintf(' The choice for the preconditioner and the solver are not compatible\n'); 0097 return 0098 end 0099 end 0100 0101 0102 npdx=nx+1; 0103 0104 [x,wx]=xwlgl(npdx); dx=derlgl(x,npdx); 0105 0106 % nov 0107 ldnov=npdx; nov=zeros(ldnov,ne); 0108 [nov]=cosnov_1d(npdx,ne,nov); 0109 noe=nov(npdx,ne); 0110 0111 % Uniform decomposition of Omega in ne elements 0112 0113 [xx,jacx,xy,ww]=mesh_1d(xa,xb,ne,npdx,nov,x,wx); 0114 0115 % ww contains the diagonal of the mass matrix 0116 0117 % G-NI Stiffness and Mass Matrices assembling 0118 if beta==0 0119 A=stiff_1d_se(npdx,ne,nov,wx,dx,jacx); 0120 A=nu*A; 0121 else 0122 A=ad_1d_se(npdx,ne,nov,nu,beta,wx,dx,jacx); 0123 end 0124 M=spdiags(ww,0,noe,noe); 0125 if gam ~=0 0126 A=A+gam*M; 0127 end 0128 % Right hand side 0129 f=(ff(xy)).*ww; 0130 0131 % Setting Dirichlet boundary conditions on the matrix on r.h.s 0132 if cb(1)=='d' 0133 f(1)=ub(1); i1=2; 0134 elseif cb(1)=='n' 0135 f(1)=f(1)-nu*ub(1);i1=1; 0136 end 0137 if cb(2)=='d'; 0138 f(noe)=ub(2);i2=noe-1; 0139 elseif cb(2)=='n' 0140 f(noe)=f(noe)+nu*ub(2);i2=noe; 0141 end 0142 0143 0144 0145 % Reduce the system to non-Dirichlet unknowns 0146 lint=(i1:i2); 0147 if cb=='dd' 0148 f(lint)=f(lint)-A(lint,1)*ub(1)-A(lint,noe)*ub(2); 0149 elseif cb=='dn' 0150 f(lint)=f(lint)-A(lint,1)*ub(1); 0151 elseif cb=='nd' 0152 f(lint)=f(lint)-A(lint,noe)*ub(2); 0153 end 0154 0155 f=f(lint); 0156 noeint=i2-i1+1; 0157 0158 if param(1)~=0 0159 % FEM Stiffness and Mass Matrices assembling 0160 [AFE,MFE,MFEd]=femp1_preco_1d(npdx,ne,nov,jacx,xy,nu,beta,gam,param); 0161 0162 % Setting Dirichlet boundary conditions on the preconditioned stiffness matrix 0163 0164 AFE=AFE(lint,lint); 0165 MFE=MFE(lint,lint); 0166 MFEd=MFEd(lint,lint); 0167 end 0168 A=A(lint,lint); 0169 M=M(lint,lint); 0170 0171 % Pre solver 0172 if param(1)==0 0173 P=speye(noeint); 0174 A1=A;f1=f; 0175 elseif param(1)==1 0176 P=AFE; 0177 A1=A;f1=f; 0178 elseif param(1)==2 0179 P=M*inv(MFE)*AFE; 0180 A1=A;f1=f; 0181 elseif param(1)==3 0182 P=M*inv(MFEd)*AFE; 0183 A1=A;f1=f; 0184 elseif param(1)==4 0185 MFE1=full(MFE);MFE1=sqrtm(MFE1);MFE1=inv(MFE1); 0186 M1=inv(sqrt(M)); 0187 P=MFE1'*AFE*MFE1; 0188 A1=M1*A*M1; 0189 f1=M1*f; 0190 elseif param(1)==5 0191 MFEd1=inv(sqrt(MFEd)); 0192 M1=inv(sqrt(M)); 0193 P=MFEd1*AFE*MFEd1; 0194 A1=M1*A*M1; 0195 f1=M1*f; 0196 end 0197 0198 % solve the linear system 0199 0200 u0=zeros(noeint,1); 0201 0202 if param(2)==1 0203 % CG to solve linear system 0204 [un,iter,res]=precg(A1, f1, u0, param(3), param(4),P,param(1)); 0205 if param(1)==4 | param(1)==5 0206 un=M1*un; 0207 end 0208 if cb=='dd' 0209 un=[ub(1);un;ub(2)]; 0210 elseif cb=='dn' 0211 un=[ub(1);un]; 0212 elseif cb=='nd' 0213 un=[un;ub(2)]; 0214 end 0215 param(11:12)=[iter,res]; 0216 0217 elseif param(2)==2 0218 % Bicgstab to solve linear system 0219 [un,iter,res]=pbcgstab_vdv(A1, f1, u0, param(3), param(4),P,param(1)); 0220 0221 if param(1)==4 | param(1)==5 0222 un=M1*un; 0223 end 0224 if cb=='dd' 0225 un=[ub(1);un;ub(2)]; 0226 elseif cb=='dn' 0227 un=[ub(1);un]; 0228 elseif cb=='nd' 0229 un=[un;ub(2)]; 0230 end 0231 param(11:12)=[iter,res]; 0232 0233 elseif param(2)==11 0234 % compute eigenvalues and the iterative condition number 0235 C=inv(P)*A1; 0236 OPTS.disp=0; 0237 d1=eigs(C,3,'LM',OPTS);d2=eigs(C,3,'SM',OPTS); 0238 kappa=max(abs(d1))/min(abs(d2)); 0239 elseif param(2)==12 0240 % compute singular values and the 2-norm condition number 0241 C=inv(P)*A1; 0242 OPTS.disp=0; 0243 d1=svds(C,3,'L',OPTS);d2=svds(C,3,0,OPTS); 0244 kappa=max(d1)/min(d2); 0245 end 0246 return 0247 0248