SHARE
TWEET

Wave equation demonstration

not_a_red_panda Sep 23rd, 2019 99 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  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
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top