• API
• FAQ
• Tools
• Archive
SHARE
TWEET

# Wave equation demonstration

not_a_red_panda Sep 23rd, 2019 203 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.
Top