Advertisement
not_a_red_panda

Wave equation demonstration

Sep 23rd, 2019
838
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 2.23 KB | None | 0 0
  1. clear
  2. clc
  3. clf
  4. %Number of Fourier coefficients
  5. N=50;
  6. %Wave velocity
  7. v=1.25;
  8. %Length
  9. L=1;
  10. %Width of pulse, 2w
  11. w = 0.2;
  12. %height of pulse
  13. h=0.4;
  14. %X axis
  15. x=linspace(0,L,100);
  16. y=zeros(1,length(x));
  17. %framerate
  18. fps=45;
  19. %Simulation time
  20. simtime=7;
  21. T=linspace(0,simtime,fps*simtime);
  22.  
  23. %Shift position of initial peak
  24. offset=0.0;
  25.  
  26. %Placeholder for Fourier and other coefficients
  27. a_n = zeros(1,N);
  28. F_n = zeros(1,N);
  29. lambda_n = zeros(1,N);
  30.  
  31. %Solution matrix
  32. soln = zeros(length(T),length(x));
  33.  
  34. %Magnitude of constant forcing, set to 0 for unforced version of problem
  35. g=-9.8;
  36.  
  37. %Damping, set to 0 for undamped version of problem
  38. b=0;
  39.  
  40. %Compute the Fourier coefficients
  41. for n=1:N
  42. bound1=0.5*L-w-offset;
  43. bound2=0.5*L-offset;
  44. bound3=0.5*L+w-offset;
  45. f1= @(s)((2/L)*sin(n*pi*s/L).*((h/w)*((s+offset)-0.5*L)+h));
  46. f2= @(s)((2/L)*sin(n*pi*s/L).*(-(h/w)*((s+offset)-0.5*L)+h));
  47. A = integral(@(s)f1(s), bound1, bound2);
  48. B = integral(@(s)f2(s), bound2, bound3);
  49. a_n(n) = A+B;
  50. F_n(n) = 2*g*(1-(-1)^n)/(n*pi);
  51. lambda_n(n) = sqrt(b^2-n^2*pi^2*v^2/L^2);
  52. end
  53.  
  54. % for n=1:N
  55. % y = y+a_n(n)*sin(n*pi*x/L);
  56. % end
  57. % plot(x,y)
  58.  
  59. for k = 1:length(T)
  60. t=T(k);
  61. y=zeros(1,length(x));
  62. for n=1:N
  63. y = y + real((1/(2*lambda_n(n)))*(a_n(n)-L^2*F_n(n)/(n^2*pi^2*v^2))*exp(-b*t)*((lambda_n(n)-b)*exp(-lambda_n(n)*t)+(b+lambda_n(n))*exp(lambda_n(n)*t))+F_n(n)*L^2/(n^2*pi^2*v^2))*sin(n*pi*x/L);
  64. end
  65. soln(k,:)=y;
  66. end
  67.  
  68. % mesh(x,T,soln)
  69. % xlabel('x');
  70. % ylabel('t');
  71.  
  72.  
  73. fig=figure;
  74. elapsed=0;
  75. nImages=length(T);
  76. for m = 1:nImages
  77. clf
  78. hold on
  79. plot(x,soln(m,:))
  80. hold off
  81. xlim([0 1])
  82. ylim([-2.25 0.5])
  83. title(['Elapsed time = ' num2str(T(m), '%.2f') ' seconds'])
  84. drawnow
  85. frame = getframe(fig);
  86. im{m} = frame2im(frame);
  87. end
  88.  
  89. close
  90. filename = 'wave2.gif';
  91. for idx = 1:nImages
  92. [A,map] = rgb2ind(im{idx},256);
  93. if idx == 1;
  94. imwrite(A,map,filename,'gif','LoopCount',Inf,'DelayTime',1.5);
  95. elseif (idx > 1) && (idx < nImages);
  96. imwrite(A,map,filename,'gif','WriteMode','append','DelayTime',0.04);
  97. elseif idx == nImages
  98. imwrite(A,map,filename,'gif','WriteMode','append','DelayTime',1.5);
  99. end
  100. end
  101. clear
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement