www.pudn.com > NumericalComputingwithMatlabCode.zip > pdegui.m, change:2003-11-30,size:12887b


function pdegui(action) 
%PDEGUI  Demonstrate solution of model partial differential equations. 
% PDEGUI demonstrates the finite difference solution of model problems 
% involving Laplace's operator: 
%      del^2_U = partial^2_U/partial_x^2 + partial^2_U/partial_y^2. 
% 
% The PDEs are: 
%      Poisson equation, del^2_U = f, f = -1. 
%      Heat equation, partial_U/dt = del^2_U - f. 
%      Wave equation, partial^2_U/partial_t^2 = del^2_U. 
%      Eigenvalue equation, del^2_U + lambda U = 0. 
% 
% The regions are:  
%      Square. 
%      L-shaped. 
%      H-shaped. 
%      Disc. 
%      Annulus. 
%      Heart. 
%      An pair of "isospectral" drums with the same eigenvalues. 
% 
% The boundary values are u = 0 outside the region.  The initial values for 
% the heat and wave equations are impulses at random points in the region. 
% You can vary the grid spacing, h.  Decreasing h increases accuracy, 
% but also increases memory requirements and computation time. 
% For the heat equation, you can set sigma = delta_t/h^2. 
% For the wave equation, you can set sigma = delta_t^2/h^2. 
% If sigma is too large, the time-stepping methods are unstable. 
% For the eigenvalue problem, you can set the eigenvalue index. 
 
if nargin == 0 
   Initialize_GUI 
   action = 'restart'; 
end 
drawnow 
 
if isequal(action,'restart') | ... 
   isequal(action,'h+') | isequal(action,'h-') 
 
   % Find the buttons 
 
   regionobj = findobj('tag','region'); 
   regionstr = get(regionobj,'string'); 
   region = regionstr{get(regionobj,'value')}; 
   eqnobj = findobj('tag','eqn'); 
   eqn = get(eqnobj,'value'); 
   hobj = findobj('tag','h'); 
   sigind = findobj('tag','sigind'); 
   indxobj = findobj('tag','indx'); 
   sigmaobj = findobj('tag','sigma'); 
   plusobj = findobj('tag','plus'); 
   minusobj = findobj('tag','minus'); 
   stopobj = findobj('string','stop'); 
 
   % Set button visibility  
 
   switch eqn 
      case {1,5}  % Poisson, Grid 
         set([sigind,indxobj,sigmaobj,plusobj,minusobj,stopobj],'vis','off'); 
      case 2      % Heat 
         set(sigind,'vis','on','string','sigma ='); 
         set(indxobj,'vis','off'); 
         set(sigmaobj,'vis','on','string','0.250','userdata',0.250); 
         set([plusobj,minusobj,stopobj],'vis','on'); 
      case 3      % Wave 
         set(sigind,'vis','on','string','sigma ='); 
         set(indxobj,'vis','off'); 
         set(sigmaobj,'vis','on','string','0.500','userdata',0.500); 
         set([plusobj,minusobj,stopobj],'vis','on'); 
      case 4      % Eigenvalue 
         set(sigind,'vis','on','string','index ='); 
         set([sigmaobj,stopobj],'vis','off'); 
         set([indxobj,plusobj,minusobj],'vis','on'); 
   end 
 
   % Adjust h 
 
   h = str2num(get(hobj,'string')); 
   if isequal(action,'h+') 
      h = 1/round(1/h+3); 
      set(hobj,'string',sprintf('1/%.0f',1/h)) 
   elseif isequal(action,'h-') & (h < 1/3) 
      h = 1/round(1/h-3); 
      set(hobj,'string',sprintf('1/%.0f',1/h)) 
   end 
 
   % Generate the mesh 
 
   [x,y] = meshgrid(-1:h:1); 
   [xv,yv,scale] = vertices(region); 
   [in,on] = inregion(x,y,xv,yv); 
   if isequal(region,'Annulus') 
      in = in - inregion(x,y,xv/2,yv/2); 
      xv = [xv NaN xv/2]; 
      yv = [yv NaN yv/2]; 
   end 
   p = find(in-on); 
   G = zeros(size(x)); 
   G(p) = 1:length(p); 
    
   % Generate the discrete Laplacian. 
 
   A = -delsq(G); 
   n = size(A,1); 
    
   % Initialize the solution 
 
   t = 0; 
   lambda = []; 
   switch eqn 
      case {2,3}      % Heat, Wave 
         p = find(G>0); 
         n = length(p); 
         q = ceil(n*rand); 
         x0 = x(p(q)); 
         y0 = y(p(q)); 
         r = (x(p)-x0).^2 + (y(p)-y0).^2; 
         if eqn == 2, c = 1; else c = 1/6; end 
         u = c*exp(-5*r); 
         v = zeros(n,1); 
      otherwise 
         u = []; 
         v = []; 
   end 
 
   % Save everything in figure's userdata. 
 
   data.eqn = eqn; data.G = G; data.A = A; data.x = x; data.y = y; 
   data.h = h; data.n = n; data.u = u; data.v = v; data.t = t; 
   data.lambda = lambda; data.xv = xv; data.yv = yv; data.scale = scale; 
   set(gcf,'userdata',data); 
 
else 
 
   % Adjust sigma or index. 
 
   eqn = get(findobj('tag','eqn'),'value'); 
   if eqn == 4 
      indxobj = findobj('tag','indx'); 
      indx = str2num(get(indxobj,'string')); 
      if isequal(action,'+') 
         indx = indx + 1; 
      elseif isequal(action,'-') 
         indx = max(1,indx - 1); 
      end 
      set(indxobj,'string',int2str(indx)); 
   else 
      sigmaobj = findobj('tag','sigma'); 
      sigma = str2num(get(sigmaobj,'string')); 
      if isequal(action,'+') 
         sigma = sigma + .001; 
      elseif isequal(action,'-') 
         sigma = sigma - .001; 
      end 
      set(sigmaobj,'string',sprintf('%6.3f',sigma)); 
   end 
 
end 
 
% Retrieve parameters 
       
data = get(gcf,'userdata'); 
A = data.A; u = data.u; v = data.v; t = data.t; scale = data.scale; 
n = data.n; h = data.h; x = data.x; y = data.y; G = data.G; 
xv = data.xv; yv = data.yv; 
 
% Begin main loop 
 
stopobj = findobj('string','stop'); 
closeobj = findobj('string','close'); 
set(stopobj,'userdata',0) 
set(closeobj,'vis','off') 
 
while get(stopobj,'userdata') == 0 
   switch eqn 
      case 1 
         % Solve sparse linear system Au = f with f = -1, 
         % scaled so that max(abs(u)) = 1. 
         f = -h^2*scale*ones(n,1); 
         u = A\f; 
         set(stopobj,'userdata',1) 
      case 2 
         % One step for heat equation 
         sigmaobj = findobj('tag','sigma'); 
         sigma = str2num(get(sigmaobj,'string')); 
         f = -h^2*scale*ones(n,1); 
         t = t + sigma*h^2; 
         u = u + sigma*(A*u - f); 
         data.t = t; 
         data.u = u; 
      case 3 
         % One step for wave equation 
         sigmaobj = findobj('tag','sigma'); 
         sigma = str2num(get(sigmaobj,'string')); 
         t = t + sqrt(sigma)*h; 
         w = u; 
         u = 2*u - v + sigma*A*u; 
         v = w; 
         data.t = t; 
         data.u = u; 
         data.v = v; 
      case 4 
         % Some eigenvalues already computed; maybe need more. 
         indxobj = findobj('tag','indx'); 
         indx = str2num(get(indxobj,'string')); 
         if indx >= n 
            indx = n-1; 
            set(indxobj,'string',num2str(n-1)) 
         end 
         while indx > length(data.lambda) 
            k = min(length(data.lambda)+10,n-1); 
            [u,lambda] = eigswatch(-(1/h^2)*A,k); 
            data.lambda = lambda; 
            data.u = u;  
         end 
         lambda = data.lambda(indx); 
         u = data.u(:,indx); 
         u = u/u(min(find(abs(u(:))==max(abs(u(:)))))); 
         set(stopobj,'userdata',1) 
      case 5 
         % Show grid, no computation 
         u = ones(n,1); 
         set(stopobj,'userdata',1) 
   end 
 
   % Map the solution onto the grid and display it. 
 
   U = zeros(size(x)); 
   p = find(G>0); 
   U(p) = u; 
   U(U>1) = 1; 
   U(U<-1) = -1; 
   set(gcf,'renderer','painters') 
   if eqn < 5 
      contourf(x,y,U,-1:1/8:1) 
   else 
      ms = max(4,min(14,ceil(100*h))); 
      plot(x(p),y(p),'.','markersize',ms,'color',[0 2/3 0]) 
   end 
   line(xv,yv,'color','black') 
   set(gca,'color',get(gcf,'color'),'xtick',[],'ytick',[]) 
   axis square 
   caxis([-1 1]) 
   switch eqn 
      case {2,3} 
         xlabel(sprintf('t = %9.4f',t)); 
      case 4 
         xlabel(sprintf('lambda(%2d) = %9.4f',indx,lambda)) 
      case 5 
         xlabel(sprintf('n = %.0f',n)); 
      otherwise 
   end 
   drawnow 
   set(gcf,'userdata',data) 
end 
set(stopobj,'vis','off') 
set(closeobj,'vis','on') 
 
 
% ------------------------------------------------------------------------ 
 
function Initialize_GUI 
clf reset 
set(gcf,'doublebuffer','on','menubar','none', ... 
   'numbertitle','off','name','PDE gui','userdata',1); 
set(gca,'pos',[.11 .11 .6 .815]) 
uicontrol( ... 
     'tag','region', ... 
     'style','list', ... 
     'units','normal', ... 
     'position',[.75 .55 .23 .35], ... 
     'string',{'Square','L','H','Disc','Annulus', ... 
        'Heart','Drum1','Drum2'}, ... 
     'fontsize',12, ... 
     'value',1, ... 
     'callback','pdegui(''restart'')') 
uicontrol( ... 
     'tag','eqn', ... 
     'style','list', ... 
     'units','normal', ... 
     'position',[.75 .30 .23 .22], ... 
     'string',{'Poisson','Heat','Wave','Eigenvalue','Grid'}, ... 
     'fontsize',12, ... 
     'value',1, ... 
     'callback','pdegui(''restart'')') 
uicontrol( ... 
     'style','text', ... 
     'units','normal', ... 
     'position',[.75 .22 .10 .06], ... 
     'background',get(gcf,'color'), ... 
     'string','   h = ', ... 
     'fontsize',12) 
uicontrol( ... 
     'tag','h', ... 
     'style','edit', ... 
     'units','normal', ... 
     'position',[.85 .22 .08 .06], ... 
     'string','1/24', ... 
     'backgroundcolor',[1 1 1], ... 
     'fontsize',12, ... 
     'callback','pdegui(''restart'')') 
uicontrol( ... 
     'tag','hminus', ... 
     'style','pushbutton', ... 
     'units','normal', ... 
     'position',[.93 .22 .03 .03], ... 
     'string','-', ... 
     'fontsize',12, ... 
     'callback','pdegui(''h-'')') 
uicontrol( ... 
     'tag','hplus', ... 
     'style','pushbutton', ... 
     'units','normal', ... 
     'position',[.93 .25 .03 .03], ... 
     'string','+', ... 
     'fontsize',12, ... 
     'callback','pdegui(''h+'')') 
uicontrol( ... 
     'tag','sigind', ... 
     'style','text', ... 
     'units','normal', ... 
     'position',[.75 .14 .10 .06], ... 
     'background',get(gcf,'color'), ... 
     'string','', ... 
     'fontsize',12) 
uicontrol( ... 
     'tag','indx', ... 
     'style','edit', ... 
     'units','normal', ... 
     'position',[.85 .14 .08 .06], ... 
     'string','1', ... 
     'userdata',1, ... 
     'backgroundcolor',[1 1 1], ... 
     'fontsize',12, ... 
     'callback','pdegui(''indx'')') 
uicontrol( ... 
     'tag','sigma', ... 
     'style','edit', ... 
     'units','normal', ... 
     'position',[.85 .14 .08 .06], ... 
     'string','0.250', ... 
     'userdata',0.250, ... 
     'backgroundcolor',[1 1 1], ... 
     'fontsize',12); 
uicontrol( ... 
     'tag','minus', ... 
     'style','pushbutton', ... 
     'units','normal', ... 
     'position',[.93 .14 .03 .03], ... 
     'string','-', ... 
     'fontsize',12, ... 
     'callback','pdegui(''-'')') 
uicontrol( ... 
     'tag','plus', ... 
     'style','pushbutton', ... 
     'units','normal', ... 
     'position',[.93 .17 .03 .03], ... 
     'string','+', ... 
     'fontsize',12, ... 
     'callback','pdegui(''+'')') 
uicontrol( ... 
     'style','pushbutton', ... 
     'units','normal', ... 
     'position',[.80 .05 .12 .06], ... 
     'string','stop', ... 
     'fontsize',12, ... 
     'userdata',0, ... 
     'callback','set(findobj(''string'',''stop''),''userdata'',1)'); 
uicontrol( ... 
     'style','toggle', ... 
     'units','normal', ... 
     'position',[.80 .05 .12 .06], ... 
     'string','close', ... 
     'fontsize',12, ... 
     'userdata',0, ... 
     'callback','close(gcf)'); 
 
 
% ------------------------------------------------------------------------ 
 
function [u,lambda] = eigswatch(A,k) 
% [u,lambda] = eigswatch(A,k) 
% Modify pointer and text while computing k smallest 
% eigenvalues and eigenvectors of sparse matrix A. 
 
ps = get(gcf,'pointer'); 
set(gcf,'pointer','watch') 
sigind = findobj('tag','sigind'); 
str = get(sigind,'string'); 
set(sigind,'string','computing') 
drawnow 
opts.disp = 0; 
[u,lambda] = eigs(A,k,0,opts); 
[lambda,p] = sort(diag(lambda)); 
u = u(:,p); 
set(sigind,'string',str) 
set(gcf,'pointer',ps) 
 
 
% ------------------------------------------------------------------------ 
 
function [xv,yv,scale] = vertices(region) 
 
% (xv,yv) = vertices of where region 
% scale = max(U) where del^2 U = -1 in region. 
 
switch region 
   case 'Square' 
      xv = [-1 1 1 -1 -1]; 
      yv = [-1 -1 1 1 -1]; 
      scale = 3.394;  
   case 'L' 
      xv = [0 1 1 -1 -1 0 0]; 
      yv = [0 0 1 1 -1 -1 0]; 
      scale = 6.702; 
   case 'H' 
      xv = [0 1 1 2 2 1 1 -1 -1 -2 -2 -1 -1 0]/2;  
      yv = [1 1 2 2 -2 -2 -1 -1 -2 -2 2 2 1 1]/2; 
      scale = 8.387; 
   case 'Disc' 
      z = exp(2*pi*i*(0:128)/128); 
      xv = real(z); 
      yv = imag(z); 
      scale = 3.960; 
   case 'Annulus' 
      z = exp(2*pi*i*(0:128)/128); 
      xv = real(z); 
      yv = imag(z); 
      scale = 30.12; 
   case 'Heart' 
      t = 0:pi/64:pi; 
      xv = sin(t); 
      yv = sin(t)-2*cos(t)/3; 
      xv = [-fliplr(xv) xv]; 
      yv = [fliplr(yv) yv]-1/3;  
      scale = 8.926; 
   case 'Drum1' 
      xv = [-3 -3 1 1 3 1 -1 -1 -3]/3; 
      yv = [-3 -1 3 1 1 -1 -1 -3 -3]/3; 
      scale = 14.96; 
   case 'Drum2' 
      xv = [-1 -3 -3 1 1 3 1 -1 -1]/3; 
      yv = [-3 -1 1 1 3 1 -1 -1 -3]/3; 
      scale = 15.35; 
end