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

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

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
'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
```