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