home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 3.ddi / DEMOS.DI$ / MEMBRANE.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  3.4 KB  |  115 lines
  1. function Lout = membrane(k,m,n,np)
  2. %MEMBRANE Generate MathWorks's logos.
  3. %
  4. %    L = MEMBRANE(k), for k <= 12, is the k-th eigenfunction of
  5. %    the L-shaped membrane.  The first three eigenfunctions have
  6. %    been shown on the covers of various MathWorks publications.
  7. %
  8. %    MEMBRANE(k), with no output parameters, plots the k-th eigenfunction.
  9. %
  10. %    MEMBRANE, with no input or output parameters, plots MEMBRANE(1).
  11. %
  12. %    L = MEMBRANE(k,m,n,np) also sets some mesh and accuracy parameters:
  13. %
  14. %      k = index of eigenfunction, default k = 1.
  15. %      m = number of points on 1/3 of boundary.  The size of
  16. %          the output is 2*m+1-by-2*m+1.  The default m = 15.
  17. %      n = number of terms in sum, default n = min(m,9).
  18. %      np = number of terms in partial sum, default np = min(n,2).
  19. %      With np = n, the eigenfunction is nearly zero on the boundary.
  20. %      With np < n, like np = 2, the boundary is not tied down.
  21. %
  22.  
  23. %    Out-of-date reference:
  24. %        Fox, Henrici & Moler, SIAM J. Numer. Anal. 4, 1967, pp. 89-102.
  25. %    Cleve Moler 4-21-85, 7-21-87, 6-30-91, 6-17-92;
  26. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  27.  
  28. if nargin < 1, k = 1; end
  29. if nargin < 2, m = 15; end
  30. if nargin < 3, n = min(m,9); end
  31. if nargin < 4, np = min(n,2); end
  32. if k > 12
  33.    error('Only the first 12 membrane eigenfunctions are available.')
  34. end
  35.  
  36. lambda = [9.6397238445, 15.19725192, 2*pi^2, 29.5214811, 31.9126360, ...
  37.     41.4745099, 44.948488, 5*pi^2, 5*pi^2, 56.709610, 65.376535, 71.057755];
  38. lambda = lambda(k);
  39. sym = [1 2 3 2 1 1 2 3 3 1 2 1];
  40. sym = sym(k);
  41.  
  42. % Part of the boundary in polar coordinates
  43. theta = (1:3*m)'/m * pi/4;
  44. r = [ones(m,1)./cos(theta(1:m)); ones(2*m,1)./sin(theta(m+1:3*m))];
  45.  
  46. % if sym = 1, alfa = [1 5 7 11 13 17 19 ... ] * 2/3, (odd, not divisible by 3)
  47. % if sym = 2, alfa = [2 4 8 10 14 16 20 ... ] * 2/3, (even, not divisible by 3)
  48. % if sym = 3, alfa = [3 6 9 12 15 18 21 ... ] * 2/3, (multiples of 3)
  49. alfa = sym;
  50. del = sym + 3*(sym == 1);
  51. for j = 2:n,
  52.    alfa(j) = alfa(j-1) + del;
  53.    del = 6-del;
  54. end;
  55. alfa = alfa * 2/3;
  56. if sym ~= 3
  57.    alf1 = (alfa(1):1:max(alfa));
  58.    alf2 = (alfa(2):1:max(alfa));
  59.    k1 = 1:4:length(alf1);
  60.    k2 = 1:4:length(alf2);
  61. else
  62.    alf1 = (alfa(1):1:max(alfa));
  63.    alf2 = [];
  64.    k1 = 1:2:length(alf1);
  65.    k2 = [];
  66. end
  67.  
  68. % Build up the matrix of fundamental solutions evaluated on the boundary.
  69. t = sqrt(lambda)*r;
  70. b1 = besselj(alf1,t);
  71. b2 = besselj(alf2,t);
  72. A = [b1(:,k1) b2(:,k2)];
  73. [ans,k] = sort([k1 k2]);
  74. A = A(:,k);
  75. A = A.*sin(theta*alfa);
  76.  
  77. % The desired coefficients are the null vector.
  78. % (lambda was chosen so that the matrix is rank deficient).  
  79. [Q,R,E] = qr(A);
  80. j = 1:n-1;
  81. c = -R(j,j)\R(j,n);
  82. c(n) = 1;
  83. c = E*c(:);
  84. if c(1) < 0, c = -c; end
  85. % The residual should be small.
  86. % res = norm(A*c,inf)/norm(abs(A)*abs(c),inf)
  87.  
  88. % Now evaluate the eigenfunction on a rectangular mesh in the interior.
  89. mm = 2*m+1;
  90. x = ones(m+1,1)*(-m:m)/m;
  91. y = (m:-1:0)'/m*ones(1,mm);
  92. r = sqrt(x.*x + y.*y);
  93. theta = atan2(y,x);
  94. theta(m+1,m+1) = 0;
  95. S = zeros(m+1,mm);
  96. r = sqrt(lambda)*r;
  97. for j = 1:np
  98.    S = S + c(j) * besselj(alfa(j),r) .* sin(alfa(j)*theta);
  99. end
  100.  
  101. L = zeros(mm,mm);
  102. L(1:m+1,:) = triu(S);
  103. if sym == 1, L = L + L' - diag(diag(L)); end
  104. if sym == 2, L = L - L'; end
  105. if sym == 3, L = L + L' - diag(diag(L)); end
  106. L = rot90(L/max(max(abs(L))),-1);
  107.  
  108. if nargout == 0
  109.    x = -1:1/m:1;
  110.    surf(x,x',L)
  111.    colormap(cool)
  112. else
  113.    Lout = L;
  114. end
  115.