LAP_1D SEM appx of the 1D b.v.p. -nu u''+gam u=f
-nu u''+gam u=f xa < x < xb
+ Dirichlet or Neumann bc
by Galerkin Numerical Integration with LGL quadrature formulas.
[xy,un,err_inf,err_l2,err_h1]=lap_1d(xa,xb,nu,gam,uex,uexx,ff,cb,ne,nx);
Input: xa, xb = extrema of computational domain Omega=(xa,xb)
nu = viscosity (constant>0)
gam = coefficient of zeroth order term (constant>=0)
uex = exact solution (uex=@(x)[uex(x)], with .*, .^, ./)
uexx = first derivative of exact solution (uexx=@(x)[uexx(x)],
with .*, .^, ./)
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
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(1) = 1: compute errors (L^inf-norm, L2-norm, H1-norm)
on the exact solution
2: no errors are computed
param(2) = 0: LG quadrature formulas with high precision degree are
used to compute norms (exact norms)
1: LGL quadrature formulas with npdx,npdy nodes are
used to compute norms (discrete norms)
(used only if param(1) == 1)
param(3) = number of nodes for high degree quadrature formula,
(used only if param(2) == 0 & param(1) == 1)
param(4) = 0: absolute errors are computed
1: relative errors are computed
(used only if param(1) == 1)
param(5) = 0: not plot the solution
1: plot the solution
param(6) = number of nodes in each element for plotting interpolated solution
(used only if param(5) ==1)
Output: xy = 1D mesh
un = numerical solution
err_inf = ||u_ex-u_N|| with respect to discrete maximum-norm
err_h1 = ||u_ex-u_N|| with respect to discrete H1-norm
err_l2 = ||u_ex-u_N|| with respect to discrete L2-norm
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 [xy,un,err_inf,err_l2,err_h1]=lap_1d(xa,xb,nu,gam,... 0002 uex,uexx,ff,cb,ne,nx,param); 0003 % LAP_1D SEM appx of the 1D b.v.p. -nu u''+gam u=f 0004 % 0005 % -nu u''+gam u=f xa < x < xb 0006 % + Dirichlet or Neumann bc 0007 % 0008 % by Galerkin Numerical Integration with LGL quadrature formulas. 0009 % 0010 % [xy,un,err_inf,err_l2,err_h1]=lap_1d(xa,xb,nu,gam,uex,uexx,ff,cb,ne,nx); 0011 % 0012 % Input: xa, xb = extrema of computational domain Omega=(xa,xb) 0013 % nu = viscosity (constant>0) 0014 % gam = coefficient of zeroth order term (constant>=0) 0015 % uex = exact solution (uex=@(x)[uex(x)], with .*, .^, ./) 0016 % uexx = first derivative of exact solution (uexx=@(x)[uexx(x)], 0017 % with .*, .^, ./) 0018 % ff = r.h.s. solution (ff=@(x)[ff(x)], with .*, .^, ./) 0019 % cb = string containing boundary conditions, for example 0020 % cb='dn' imposes Dirichlet in xa and Neumann in xb 0021 % ne = number of elements (equally spaced) 0022 % nx = polynomial degree in each element (the same in each element) 0023 % to be set only if p=4, otherwise, nx=p; 0024 % param(1) = 1: compute errors (L^inf-norm, L2-norm, H1-norm) 0025 % on the exact solution 0026 % 2: no errors are computed 0027 % param(2) = 0: LG quadrature formulas with high precision degree are 0028 % used to compute norms (exact norms) 0029 % 1: LGL quadrature formulas with npdx,npdy nodes are 0030 % used to compute norms (discrete norms) 0031 % (used only if param(1) == 1) 0032 % param(3) = number of nodes for high degree quadrature formula, 0033 % (used only if param(2) == 0 & param(1) == 1) 0034 % param(4) = 0: absolute errors are computed 0035 % 1: relative errors are computed 0036 % (used only if param(1) == 1) 0037 % param(5) = 0: not plot the solution 0038 % 1: plot the solution 0039 % param(6) = number of nodes in each element for plotting interpolated solution 0040 % (used only if param(5) ==1) 0041 % 0042 % Output: xy = 1D mesh 0043 % un = numerical solution 0044 % err_inf = ||u_ex-u_N|| with respect to discrete maximum-norm 0045 % err_h1 = ||u_ex-u_N|| with respect to discrete H1-norm 0046 % err_l2 = ||u_ex-u_N|| with respect to discrete L2-norm 0047 % 0048 % Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0049 % "Spectral Methods. Fundamentals in Single Domains" 0050 % Springer Verlag, Berlin Heidelberg New York, 2006. 0051 0052 % Written by Paola Gervasio 0053 % $Date: 2007/04/01$ 0054 0055 xy=0;un=0;err_inf=0;err_l2=0;err_h1=0; 0056 0057 0058 % 0059 npdx=nx+1; 0060 0061 [x,wx]=xwlgl(npdx); dx=derlgl(x,npdx); 0062 0063 % nov 0064 ldnov=npdx; nov=zeros(ldnov,ne); 0065 [nov]=cosnov_1d(npdx,ne,nov); 0066 noe=nov(npdx,ne); 0067 0068 % Uniform decomposition of Omega in ne elements 0069 0070 [xx,jacx,xy,ww]=mesh_1d(xa,xb,ne,npdx,nov,x,wx); 0071 0072 % ww contains the diagonal of the mass matrix 0073 0074 % Stiffness Matrix assembling 0075 A=stiff_1d_se(npdx,ne,nov,wx,dx,jacx); A=A*nu+gam*diag(ww); 0076 0077 % Setting Dirichlet boundary conditions on the matrix 0078 0079 if cb(1)=='d'; A(1,:)=zeros(1,noe); A(1,1)=1; end 0080 if cb(2)=='d'; A(noe,:)=zeros(1,noe); A(noe,noe)=1; end 0081 0082 % Right hand side 0083 f=zeros(noe,1); 0084 f=ff(xy).*ww; 0085 0086 % Setting Dirichlet boundary conditions on both matrix and r.h.s 0087 0088 if cb(1)=='d'; 0089 f(1)=uex(xy(1)); 0090 else 0091 f(1)=f(1)-nu*uexx(xy(1)); 0092 end 0093 if cb(2)=='d'; 0094 f(noe)=uex(xy(noe)); 0095 else 0096 f(noe)=f(noe)+nu*uexx(xy(noe)); 0097 end 0098 0099 % direct solution of the linear system 0100 0101 un=A\f; 0102 0103 if(param(1)==1) 0104 % Evaluate exact solution to compute errors 0105 nq=param(3); fdq=param(2); err_type=param(4); 0106 0107 u=uex(xy); 0108 0109 % 0110 0111 err=u-un; 0112 0113 % ||u_ex-u_N|| with respect to discrete maximum-norm 0114 0115 if err_type==0 0116 err_inf=norm(err,inf); 0117 else 0118 err_inf=norm(err,inf)/norm(u,inf); 0119 end 0120 0121 0122 % ||u_ex-u_N|| with respect to H1 norm 0123 % fdq=0 LG quadrature formula with nq nodes in each element 0124 % fdq=1 LGL quadrature formula with npdx nodes in each element 0125 [err_h1]=normah1_1d(fdq, nq, err_type, un, uex, uexx,... 0126 x, wx, dx, xx, jacx, xy,nov); 0127 0128 % ||u_ex-u_N|| with respect to L2 norm 0129 0130 [err_l2]=normal2_1d(fdq, nq, err_type, un, uex, ... 0131 x, wx, xx, jacx, xy,nov); 0132 end 0133 0134 if(param(5)==1) 0135 % interpolate numerical solution to give a "good" plot 0136 0137 fig=figure(..., 0138 'Name','SEM-NI solution of -nu u'''' +gam u=f xa < x < xb',... 0139 'Visible','on'); 0140 ha=plot_sem_1d(fig,ne,x,wx,xx,jacx,xy,nov,un,param(6)); 0141 end 0142 0143 0144 return