ell_1d

 ELL_1D   Numerical solution of the 1D elliptic boundary value problem (4.3.17) pag. 208, CHQZ2

  -nu d^2 u/dx^2 +beta du/dx+ gam u=1  in (-1,1)

   + Dirichlet or Neumann boundary conditions

 by Galerkin Numerical Integration with LGL quadrature formulas.

  [xy,un]=ell_1d(xa,xb,nu,beta,gam,ff,ub,ne,nx);

 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)
      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)

 Output: xy = 1D mesh
         un = numerical solution 


 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]=ell_1d(xa,xb,nu,beta,gam,ff,cb,ub,ne,nx);
0002 % ELL_1D   Numerical solution of the 1D elliptic boundary value problem (4.3.17) pag. 208, CHQZ2
0003 %
0004 %  -nu d^2 u/dx^2 +beta du/dx+ gam u=1  in (-1,1)
0005 %
0006 %   + Dirichlet or Neumann boundary conditions
0007 %
0008 % by Galerkin Numerical Integration with LGL quadrature formulas.
0009 %
0010 %  [xy,un]=ell_1d(xa,xb,nu,beta,gam,ff,ub,ne,nx);
0011 %
0012 % Input: xa, xb = extrema of computational domain Omega=(xa,xb)
0013 %      nu   = viscosity (constant>0)
0014 %      beta = coefficient of first order term (constant)
0015 %      gam = coefficient of zeroth order term (constant>0)
0016 %      cb = string containing boundary conditions, for example
0017 %           cb='dn' imposes Dirichlet in xa and Neumann in xb,
0018 %           possible b.c.: Dirichlet, Neumann
0019 %      ub = 2 components array with Dirichlet, Neumann (pure derivative)
0020 %           data in xa and in xb
0021 %      ne = number of elements (equally spaced)
0022 %      nx = polynomial degree in each element (the same in each element)
0023 %
0024 % Output: xy = 1D mesh
0025 %         un = numerical solution
0026 %
0027 %
0028 % Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang,
0029 %                    "Spectral Methods. Fundamentals in Single Domains"
0030 %                    Springer Verlag, Berlin Heidelberg New York, 2006.
0031 
0032 %   Written by Paola Gervasio
0033 %   $Date: 2007/04/01$
0034 
0035 
0036 %
0037 npdx=nx+1; 
0038 
0039 [x,wx]=xwlgl(npdx); dx=derlgl(x,npdx);
0040 
0041 % nov
0042 ldnov=npdx; nov=zeros(ldnov,ne);
0043 [nov]=cosnov_1d(npdx,ne,nov);
0044 noe=nov(npdx,ne);
0045 
0046 % Uniform decomposition of  Omega in ne elements
0047 
0048 [xx,jacx,xy,ww]=mesh_1d(xa,xb,ne,npdx,nov,x,wx);
0049 
0050 % Variable coefficients evaluation
0051 
0052 
0053 % Matrix assembling
0054 A=ell_1d_se(npdx,ne,nov,nu,beta,gam,wx,dx,jacx);
0055 
0056 % Right hand side
0057 f=ff(xy).*ww;
0058 
0059 
0060 % Setting Dirichlet boundary conditions on the matrix
0061 
0062 if cb(1)=='d'
0063 f(1)=ub(1); i1=2;
0064 elseif cb(1)=='n'
0065 f(1)=f(1)-nu*ub(1);i1=1;
0066 end
0067 if cb(2)=='d';
0068 f(noe)=ub(2);i2=noe-1;
0069 elseif cb(1)=='n'
0070 f(noe)=f(noe)+nu*ub(2);i2=noe;
0071 end
0072 
0073 % Reduce the system to non-Dirichlet unknowns
0074 lint=(i1:i2);
0075 if cb=='dd'
0076 f(lint)=f(lint)-A(lint,1)*ub(1)-A(lint,noe)*ub(2);
0077 elseif cb=='dn'
0078 f(lint)=f(lint)-A(lint,1)*ub(1);
0079 elseif cb=='nd'
0080 f(lint)=f(lint)-A(lint,noe)*ub(2);
0081 end
0082 
0083 A=A(lint,lint);
0084 f=f(lint);
0085 noeint=i2-i1+1;
0086 
0087 % MEG on the linear system
0088 
0089 un=A\f;
0090 
0091 if cb=='dd'
0092 un=[ub(1);un;ub(2)];
0093 elseif cb=='dn'
0094 un=[ub(1);un];
0095 elseif cb=='nd'
0096 un=[un;ub(2)];
0097 end
0098 
0099 return