home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 3.ddi / DEMOS.DI$ / DELSQSHO.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  1.6 KB  |  65 lines
  1. function delsqshow(n,Rlist)
  2. %DELSQSHOW    Show discrete Laplacian on various grids.
  3. %    delsqshow(n,Rlist)
  4. %    Rectangular grid is n-by-n.
  5. %    Rlist is a string of letters denoting regions, taken
  6. %    from the list 'SNLCDAHB'.  See numgrid.
  7.  
  8. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  9.  
  10. clf
  11.  
  12. % Grid size
  13. if nargin < 1, n = 32; end
  14. if nargin < 2, Rlist = 'SNLCDAHB'; end
  15.    
  16. % Loop over regions known to "numgrid"
  17. for R = Rlist
  18.  
  19.     % Generate and display the grid.
  20.     G = numgrid(R,n);
  21.     spy(G)
  22.     title('A finite difference grid')
  23.     g = numgrid(R,12)
  24.     disp('pause'), disp(' '), pause
  25.  
  26.     % Generate and display the discrete Laplacian.
  27.     D = delsq(G);
  28.     spy(D)
  29.     title('The 5-point Laplacian')
  30.     disp('pause'), disp(' '), pause
  31.  
  32.     % Number of interior points
  33.     N = sum(G(:)>0);
  34.  
  35.     % Solve the Dirichlet boundary value problem:
  36.     %    delsq(u) = 1 in the interior,
  37.     %    u = 0 on the boundary.
  38.     disp('Solving the sparse linear system ...')
  39.     rhs = ones(N,1);
  40.     tic
  41.     if R == 'N'
  42.         % For nested dissection, turn off minimum degree ordering.
  43.         spparms('autommd',0)
  44.         u = D\rhs;
  45.         spparms('autommd',1)
  46.     else
  47.         u = D\rhs;
  48.     end
  49.     toc
  50.  
  51.     % Map the solution onto the grid and display it.
  52.     U = G;
  53.     U(G>0) = full(u(G(G>0)));
  54.     clabel(contour(U));
  55.     prism
  56.     axis('square'), axis('ij')
  57.     disp('pause'), disp(' '), pause
  58.  
  59.     colormap((cool+1)/2);
  60.     mesh(U)
  61.     axis([0 n 0 n 0 max(max(U))])
  62.     axis('square'), axis('ij')
  63.     disp('pause'), disp(' '), pause
  64. end
  65.