home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l460 / 3.ddi / DEMOS.DI$ / VIBES.M < prev    next >
Encoding:
Text File  |  1993-03-07  |  2.8 KB  |  105 lines
  1. %VIBES    Movie of the vibrating L-shaped membrane.
  2. echo off
  3. %    C. B. Moler, 6-30-91, 8-27-92.
  4. %    Copyright (c) 1984-93 by The MathWorks, Inc.
  5.  
  6. clf
  7. drawnow
  8. echo on
  9. clc
  10.  
  11. % This demonstration solves the wave equation for the vibrations
  12. % of an L-shaped membrane.  The solution is expressed as a linear
  13. % combination, with time-dependent coefficients, of two-dimensional
  14. % spatial eigenfunctions.  The eigenfunctions have been pre-computed,
  15. % using the function MEMBRANE, and saved in a file.  The first of
  16. % these eigenfunctions, the fundamental mode, is the MathWorks logo.
  17. %
  18. % The L-shaped geometry is of particular interest mathematically because
  19. % the stresses approach infinity near the reentrant corner.  Conventional
  20. % finite difference and finite element methods require considerable
  21. % time and storage to achieve reasonable accuracy.  The approach used
  22. % here employs Bessel functions with fractional order to match the
  23. % corner singularity.
  24. %
  25. % As the solution is computed and displayed at each of 12 time steps,
  26. % snapshots of the graphics window are saved in a large matrix.
  27. % The total storage required is proportional to the area of the graphics
  28. % window.  About 2 megabytes is required for the default window.
  29.  
  30. echo off
  31.  
  32. % Load eigenfunction data.
  33. load vibesdat
  34. n = max(size(L1));
  35. nh = fix(n/2);
  36. x = (-nh:nh)/nh;
  37. clf
  38. contour(x,x,L1,18), prism
  39. disp('Press any key to proceed ...')
  40. pause
  41.  
  42. % Get coefficients from eigenfunctions.
  43. clear c
  44. for k = 1:12
  45.    eval(['c(k) = L' num2str(k) '(24,13)/3;'])
  46. end
  47.  
  48. % Set graphics parameters.
  49. clf
  50. axis([-1 1 -1 1 -1 1]);
  51. caxis(26.9*[-1.5 1]);
  52. colormap(hot);
  53. format compact
  54. hold on
  55.  
  56. % Generate the movie.
  57. delt = 0.1;
  58. nframes = 12;
  59. M = moviein(nframes);
  60. for k = 1:nframes
  61.    disp(k)
  62.  
  63.    % Coefficients
  64.    t = k*delt;
  65.    s = c.*sin(sqrt(lambda)*t);
  66.  
  67.    % Amplitude
  68.    L = s(1)*L1 + s(2)*L2 + s(3)*L3 + s(4)*L4 + s(5)*L5 + s(6)*L6 + ...
  69.        s(7)*L7 + s(8)*L8 + s(9)*L9 + s(10)*L10 + s(11)*L11 + s(12)*L12;
  70.  
  71.    % Velocity
  72.    s = s .* lambda;
  73.    V = s(1)*L1 + s(2)*L2 + s(3)*L3 + s(4)*L4 + s(5)*L5 + s(6)*L6 + ...
  74.        s(7)*L7 + s(8)*L8 + s(9)*L9 + s(10)*L10 + s(11)*L11 + s(12)*L12;
  75.  
  76.    % Surface plot of height, colored by velocity.
  77.    % Cut out the reentrant corner
  78.    V(1:nh,1:nh) = NaN*ones(nh,nh);
  79.    cla
  80.    surf(x,x,L,V);
  81.  
  82.    % Capture the frame
  83.    M(:,k) = getframe;
  84. end
  85. hold off
  86.  
  87. echo on
  88.  
  89. % Now the movie can be played back at roughly 12 frames per second.
  90.  
  91. movie(M,20,12);
  92.  
  93. echo off
  94.  
  95. % Push buttons
  96.  
  97. callback = [ ...
  98.    'set(uip,''back'',[1 2/3 1/3]);' ...
  99.    'movie(M,20,12); ' ...
  100.    'set(uip,''back'',''default'');' ];
  101. uip = uicontrol('style','push','pos',[20 10 60 40],'string','Replay', ...
  102.    'call',callback);
  103. uiq = uicontrol('style','push','pos',[100 10 60 40],'string','Done', ...
  104.    'call','clf; ');
  105.