STAT_SCALAR_HYP Numerical solution of a stationary scalar linear hyperbolic equation,
formula (3.7.1a), pag. 145, CHQZ2
beta du/dx=f(x) in (xa,xb)
inflow cond.: u(xa)=ul
Space discretization: 1. Galerkin-Numerical Integration with LGL formulas,
2. Legendre collocation (strong form),
3. penalty,
[x,u,err]=stat_scalar_hyp(xa,xb,beta,f,uex,ul,nx,param);
Input: xa,xb = extrema of space domain Omega=(xa,xb)
beta = transport coefficient
uex = exact solution (uex=@(x)[uex(x)], with .*, .^, ./)
ul = inflow condition (scalar datum)
nx = spectral polynomial degree, the number of LGL points is
nx+1
param(1)= space disc. scheme: 1=strong (3.7.2), pag. 146 CHQZ2,
2=GNI+ weak b.c. (3.7.5), pag. 147
3=penalty (3.7.2a)+(3.7.9), pag. 148
param(2)=tau : penalty parameter (used only if param(1)==3)
Output: u = solution
err= L_\infty norm error between numerical and exact solution
Possible input data:
xa=-1;xb=1;
beta=1.5;
param(1)=2;param(2)=0.75
nx=8;
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 [x,u,err]=stat_scalar_hyp(xa,xb,beta,f,uex,ul,nx,param) 0002 % STAT_SCALAR_HYP Numerical solution of a stationary scalar linear hyperbolic equation, 0003 % formula (3.7.1a), pag. 145, CHQZ2 0004 % 0005 % beta du/dx=f(x) in (xa,xb) 0006 % inflow cond.: u(xa)=ul 0007 % 0008 % Space discretization: 1. Galerkin-Numerical Integration with LGL formulas, 0009 % 2. Legendre collocation (strong form), 0010 % 3. penalty, 0011 % 0012 % 0013 % [x,u,err]=stat_scalar_hyp(xa,xb,beta,f,uex,ul,nx,param); 0014 % 0015 % Input: xa,xb = extrema of space domain Omega=(xa,xb) 0016 % beta = transport coefficient 0017 % uex = exact solution (uex=@(x)[uex(x)], with .*, .^, ./) 0018 % ul = inflow condition (scalar datum) 0019 % nx = spectral polynomial degree, the number of LGL points is 0020 % nx+1 0021 % param(1)= space disc. scheme: 1=strong (3.7.2), pag. 146 CHQZ2, 0022 % 2=GNI+ weak b.c. (3.7.5), pag. 147 0023 % 3=penalty (3.7.2a)+(3.7.9), pag. 148 0024 % param(2)=tau : penalty parameter (used only if param(1)==3) 0025 % 0026 % Output: u = solution 0027 % err= L_\infty norm error between numerical and exact solution 0028 % 0029 % Possible input data: 0030 % xa=-1;xb=1; 0031 % beta=1.5; 0032 % param(1)=2;param(2)=0.75 0033 % nx=8; 0034 % 0035 % Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0036 % "Spectral Methods. Fundamentals in Single Domains" 0037 % Springer Verlag, Berlin Heidelberg New York, 2006. 0038 0039 % Written by Paola Gervasio 0040 % $Date: 2007/04/01$ 0041 0042 0043 0044 npdx=nx+1; 0045 [x,wx]=xwlgl(npdx); 0046 [dx]=derlgl(x,npdx); 0047 jac=(xb-xa)*5.d-1; 0048 x=x*jac+(xa+xb)*5.d-1; 0049 0050 0051 tau=param(2)/(jac*wx(1)); 0052 0053 % Construction of matrix A 0054 0055 if param(1)==1 0056 % strong collocation (3.7.2) 0057 A=beta*diag(wx)*dx; 0058 A(1,:)=zeros(1,npdx); A(1,1)=1; 0059 0060 0061 elseif param(1)==2 0062 % GNI+weak b.c. (3.7.5) 0063 A=-beta*dx'*diag(wx); 0064 A(npdx,npdx)=A(npdx,npdx)+beta; 0065 0066 elseif param(1)==3 0067 % penalty (3.7.2a) +(3.7.9) 0068 A=beta*diag(wx)*dx; 0069 A(1,1)=A(1,1)+tau*beta*wx(1)*jac; 0070 0071 end 0072 0073 [L,U,P]=lu(A); 0074 0075 %right hand side 0076 b=f(x).*wx*jac; 0077 0078 0079 % inflow condition 0080 if param(1)==1 0081 b(1)=ul; 0082 elseif param(1)==2 0083 b(1)=b(1)+beta*ul; 0084 elseif param(1)==3 0085 b(1)=b(1)+tau*beta*ul*wx(1)*jac; 0086 end 0087 0088 % linear system solution 0089 u=U\(L\(P*b)); 0090 0091 0092 % error 0093 0094 ue=uex(x); 0095 err=norm(u-ue); 0096 0097 return