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. OK, I Understand
 
Top