• API
• FAQ
• Tools
• Archive
SHARE
TWEET

Untitled

a guest Jun 18th, 2019 74 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. clear all; close all;
2.
3. % --- Define constants and initial condition
4. alpha1 = 0.065;
5. alpha2 = 0.012;
6. L = 50.0; % length of domain in x direction in cm
7. tmax = 1000.0;  % end time in seconds
8.
9.
10.
11. tNodes = 10000.0;
12. xNodes = 200.0;
13. %%
14. dx = 1.0/xNodes;
15. dt = 1.0/tNodes;
16. %%
17. r1 = (alpha1*dt)/(dx^2);
18. r2 = 1 - 2*r1;
19.
20. rStar = (alpha2*dt)/(dx^2);
21. r2Star = 1 - 2*rStar;
22. %%
23. sigma = 10; %FWHM in cm
24. % --- Loop over time steps
25. tArray = linspace(0,tmax,tNodes);
26. xArray = linspace(-L,L,xNodes);
27.
28.
29. f = @(x) (1/sqrt(2*pi*sigma.^2))*exp(-x.^2/(4*sigma.^2));
30. u = f(xArray);  % initial condition
31.
32. %%%%%%%%%%%%%%%%%%%
33. xSize = length(xArray);
34. tSize = length(tArray);
35. uArray = zeros(tSize, xSize);
36.
37.
38. %Initializing u Array
39. uArray(1,:) = u(:);
40.
41. %Creating laplacian matrix
42. %first, creating r-array for varying coefficient
43. rArray = ones(1,200);
44. rArray(1:89)    = r2;
45. rArray(90:94)  = r2Star;
46. rArray(95:105)  = r2;
47. rArray(106:110)  = r2Star;
48. rArray(111:200) = r2;
49.
50. rArray = transpose(rArray);
51.
52. rArray2 = ones(1,200);
53. rArray2(1:89)    = r1;
54. rArray2(90:94)  = rStar;
55. rArray2(95:105)  = r1;
56. rArray2(106:110)  = rStar;
57. rArray2(111:200) = r1;
58.
59. diagR = eye(200);
60. i = 1;
61. while i <= 200
62.     diagR(i,i) = rArray(i);
63.     i = i + 1;
64. end
65.
66. Atemp = diagR;
67.
68. ASupDiag = zeros(200,200);
69. ASubDiag = zeros(200,200);
70.
71.
72. j = 1;
73. while j <= 100
74.     ASupDiag(j, j+1) = rArray2(j);
75.     ASubDiag(j + 1,j) = rArray2(j);
76.     j = j + 1;
77. end
78.
79. % ASubDiag(201,:) = [];
80. % ASupDiag(:,201) = [];
81.
82. A = Atemp + ASubDiag + ASupDiag;
83. A(1,1) = 1;
84. A(end,end) = 1;
85. A(end, end - 1) = 0;
86. A(1,2) = 0;
87.
88. Prop = A;
89. %%%% Finished creating propagator
90.
91.
92. %%
93.
94. for m=1:tNodes - 1
95.     uTemp = reshape(uArray(m,:),[],1);
96.
97.     uArray(m + 1, :) = A*uTemp;
98. end
99.
100.
101. figure
102. colormap hsv
103. surf(xArray, tArray, uArray,'FaceColor','interp',...
104.    'EdgeColor','none',...
105.    'FaceLighting','gouraud')
106. ylim([0 100])
107. xlabel('x (cm)')
108. ylabel('t (seconds)')
109. zlabel('$\rho$', 'interpreter', 'latex')
110. camlight left
111.
112.
113. %% Integration
114.
115. NormVec = uArray(1,:);
116. norm = sum(NormVec);
117. uArray = uArray/norm;
118.
119.
120. for j = 1:64000
121.     sumArray = reshape(uArray(j, 45:55), [], 1);
122.     sum1 = sum(sumArray);
123.     if sum1 < 0.01
124.         timestamp = tArray(j);
125.         break
126.     else
127.         continue
128.     end
129. end
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