www.pudn.com > NumericalComputingwithMatlabCode.zip > membranetx.m, change:2003-12-25,size:6495b


function [L,lambda] = membranetx(k,m,n,np) 
%MEMBRANETX  Textbook version of MEMBRANE, eigenfunctions of L-membrane. 
% 
%   L = MEMBRANETX(k) is the k-th eigenfunction of the L-shaped membrane. 
%   [L,lambda] = MEMBRANETX(k) also returns the k-th eigenvalue. 
% 
%   L = MEMBRANETX(k,m,n,np) sets some mesh and accuracy parameters: 
% 
%     k = index of eigenfunction, default k = 1. 
%     m = number of points on one edge of one square. 
%         The output L is 2*m+1-by-2*m+1.  The default m = 30. 
%     n = number of terms in sum, default n = min(m,20). 
%     np = number of terms in partial sum, default np = n. 
%     With np = n, the eigenfunction is zero on the boundary. 
%     With np < n, such as np = 2, the boundary is not tied down. 
% 
%   L = ROT90(MEMBRANETX(1,15,9,2),-1) is the MathWorks logo. 
 
% Default parameters 
 
if nargin < 1, k = 1; end 
if nargin < 2, m = 30; end 
if nargin < 3, n = min(m,20); end 
if nargin < 4, np = n; end 
 
% Compute eigenvalue and symmetry class. 
% sym = 1, symmetric about center line 
% sym = 2, antisymmetric about center line 
% sym >= 3, eigenvalue of the square, reflected into other squares 
 
[lambda,sym] = membraneval(k,m,n); 
if m == 1, L = lambda(k); return, end 
 
% The null vector from the SVD of the boundary matrix gives coefficients. 
 
[sigma,c,alfa] = membranesvd(lambda,sym,m,n); 
 
% Evaluate the eigenfunction on a square grid. 
 
L = membranefun(lambda,sym,c,alfa,m,n,np); 
 
% ------------------------------ 
 
function [lambda,sym] = membraneval(k,m,n); 
% MEMBRANEVAL 
% [lambda,sym] = membraneval(k,m,n) is the k-th eigenvalue of 
% the L shaped membrane, and its symmetry class. 
% m = number of points on one edge of one square. 
% n = number of terms in sum. 
 
persistent lambdas syms 
 
if isempty(lambdas) & exist('membrane.mat') 
   % Load precomputed eigenvalues 
   load membrane.mat 
end 
 
if length(lambdas) < k 
   % Compute eigenvalues beyond those already computed. 
   % Algorithm: 
   % Use direct search to get near local minima of membranesvd(lambda). 
   % Then use "fmintx" to home in on the minimizers. 
   % The step size delta controls the direct search. 
   % Increasing delta decreases computer time, but might miss some eigenvalues.  
   % kmax = number of eigenvalues. 
   % delta = search increment. 
   % tol = tolerance for fmintx 
   kmax = k; 
   delta = .01; 
   tol = 1.e-12; 
   k = length(lambdas); 
   if k == 0 
      lambdas(1) = fmintx(@membranesvd,9.6,9.7,tol,1,m,n); 
      syms(1) = 1; 
      k = 1; 
      fprintf(1,'%4.0d %18.12f %4.0d\n',k,lambdas(k),syms(k)) 
   end 
   xstart = delta*floor(lambdas(k)/delta); 
   x = [0 0 xstart]; 
   f = zeros(3,3); 
 
   % Look for x so that f(x) < both f(x-delta) and f(x+delta). 
   while k < kmax 
      x(1:2) = x(2:3); 
      x(3) =  x(3) + delta; 
      for s = 1:3    % Symmetry class. 
         f(s,1:2) = f(s,2:3); 
         f(s,3) = membranesvd(x(3),s,m,n); 
         if f(s,2) < f(s,1) & f(s,2) < f(s,3); 
            lam = fmintx(@membranesvd,x(1),x(3),tol,s,m,n); 
            if s < 3 
               mult = 1; 
            else 
               % Multiple eigenvalues are integer multiples of pi^2 
               p = round(lam/pi^2); 
               lam = p*pi^2; 
               [i,j] = ndgrid(1:sqrt(p)); 
               mult = sum(p == i(:).^2+j(:).^2); 
            end 
            for mu = 1:mult 
               k = k+1; 
               lambdas(k,1) = lam; 
               syms(k,1) = s+mu-1; 
               fprintf(1,'%4.0d %18.12f %4.0d\n',k,lambdas(k),syms(k)) 
               pause(0) 
            end 
         end 
      end 
   end 
end 
 
[lambdas,p] = sort(lambdas); 
syms = syms(p); 
% save membrane lambdas syms 
lambda = lambdas(k); 
sym = syms(k); 
 
% ------------------------------ 
 
function [sigma,c,alfa] = membranesvd(lambda,sym,m,n) 
% MEMBRANESVD 
% Evaluate fundamental solutions on boundary of L-shaped region. 
% sigma = membranesvd(lambda,s,m,n) is the smallest singular value of the 
% matrix obtained by evaluating n fundamental solutions with symmetry class 
% s at 3*m+1 points on the boundary of the L.  If lambda is chosen to give 
% a local minima of this function, the resulting null vector, c, provides 
% coeffients for a linear combination over the entire region that nearly 
% vanishes on the boundary. 
% 
% Input: 
%   lambda = eigenvalue parameter, vary this to minimize the resulting sigma. 
%   sym = symmetry class. 
%   m = number of points on edge of one square. 
%   n = number of fundamental solutions. 
% 
% Output: 
%   sigma = smallest singular value. 
%   c = null vector = coefficients. 
%   alfa = 1-by-n vector of Bessel function orders for given symmetry. 
 
% Bessel function orders. 
% sym = 1, alfa = (2/3) * [1 5 7 11 13 ... ], (odd, not divisible by 3) 
% sym = 2, alfa = (2/3) * [2 4 8 10 14 ... ], (even, not divisible by 3) 
% sym >= 3, alfa = [2 4 6 8 10 ... ] = even integers 
 
switch sym 
   case {1,2} 
      j = (sym:2:3*n); 
      j(mod(j,3)==0) = []; 
      alfa = (2/3)*j; 
   otherwise 
      alfa = 2*(1:n); 
end 
 
% Use polar coordinates to describe three-eighths of the boundary. 
 
x = [ones(m,1); (m:-1:-m)'/m]; 
y = [(0:m-1)'/m; ones(2*m+1,1)]; 
theta = atan2(y,x); 
r = sqrt(x.^2 + y.^2); 
 
% Evaluate the fundamental solutions on the boundary. 
% A is a (3*m+1)-by-n matrix. 
 
A = besselj(alfa,sqrt(lambda)*r).*sin(theta*alfa); 
 
% Scale to make columns comparable. 
 
scale = diag(sparse(1./sqrt(sum(A.*A)))); 
A = A*scale; 
 
% Compute SVD and obtain coefficients from null vector(s). 
 
[U,S,V] = svd(A,0); 
if sym > 3, n = n-(sym-3); end    % Multiple eigenvalue 
sigma = S(n,n); 
c = scale*V(:,n); 
 
% ------------------------------ 
 
function L = membranefun(lambda,sym,c,alfa,m,n,np) 
% MEMBRANEFUN  Evaluate the eigenfunction on a square grid. 
% L = membranefun(lambda,sym,c,alfa,m,n,np) 
% Used by MEMBRANETX. 
 
[x,y] = meshgrid((-m:m)/m,(m:-1:0)'/m); 
r = sqrt(x.*x + y.*y); 
theta = atan2(y,x); 
theta(m+1,m+1) = 0; 
S = zeros(m+1,2*m+1); 
for j = 1:np 
   S = S + c(j)*besselj(alfa(j),sqrt(lambda)*r).*sin(alfa(j)*theta); 
end 
S = S/S(min(find(abs(S(:)) == max(abs(S(:)))))); 
L = zeros(2*m+1,2*m+1); 
switch sym 
   case 1 
      L(1:m+1,:) = triu(S); 
      L = L + L' - diag(diag(L)); 
   case 2 
      L(1:m+1,:) = triu(S); 
      L = L - L'; 
   otherwise 
      L(1:m,1:m) = S(1:m,1:m); 
      L(m+2:2*m+1,1:m) = -flipud(L(1:m,1:m)); 
      L(1:m,m+2:2*m+1) = -fliplr(L(1:m,1:m)); 
end