Home > Src > Elliptic_1d > lap_1d.m

lap_1d

PURPOSE ^

LAP_1D SEM appx of the 1D b.v.p. -nu u''+gam u=f

SYNOPSIS ^

function [xy,un,err_inf,err_l2,err_h1]=lap_1d(xa,xb,nu,gam,uex,uexx,ff,cb,ne,nx,param);

DESCRIPTION ^

 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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Fri 21-Sep-2007 10:07:00 by m2html © 2003