Advertisement
Guest User

Untitled

a guest
May 21st, 2018
124
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 1.98 KB | None | 0 0
  1. %
  2. % This script uses Godunov scheme to
  3. % simulate the traffic flow problem
  4. %
  5.  
  6. clear all
  7.  
  8. global rho_max;
  9. global u_max;
  10.  
  11. %
  12. % problem parameters
  13. %
  14. rho_max=200.0
  15. u_max=60.0
  16. rho_L=20.0
  17.  
  18. %
  19. % we will consider x in range [-2,2], and t in range [0,2]
  20. %
  21. T=4/60;
  22. XL=-2;
  23. XR=2;
  24.  
  25. %
  26. % discretization
  27. %
  28. delta_x = 4/400
  29. delta_t = 0.8*delta_x/u_max
  30.  
  31. delt_by_delx = delta_t/delta_x;
  32.  
  33. xs = XL:delta_x:XR;
  34. ts = 0:delta_t:T;
  35. Nx = length(xs);
  36. Nt = length(ts);
  37. rho = zeros(Nt,Nx);
  38.  
  39. %
  40. % put your initial conditions here
  41. %
  42. i = 1;
  43. Nx0 = -XL/delta_x;
  44. rho(1,:) = 100;
  45. %rho(i,:) = [rho_max*ones(1, Nx0), rho_L*ones(1, Nx-Nx0)];
  46. %
  47. % a few sample trajectories
  48. %
  49.  
  50. traj(1,:) = [-2.0 -1.0 1.0];
  51. Ntraj = length(traj);
  52.  
  53. for n=2:Nt
  54.  
  55. % incoming flux at x = XL
  56. flux_minus = numflux(rho(n-1,1), rho(n-1,1));
  57.  
  58. for i=1:Nx
  59.  
  60. % compute discrete fluxes F(xi-1/2) and F(xi+1/2)
  61. % assume the disturbance never reaches the domain boundaries
  62. % Note: F(xi-1/2) is taken from the previous step of the loop
  63.  
  64. rho_i = rho(n-1, i);
  65.  
  66. if (i==Nx)
  67. rho_i_plus_1 = rho_i;
  68. else
  69. rho_i_plus_1 = rho(n-1,i+1);
  70. end
  71.  
  72. flux_plus = numflux(rho_i, rho_i_plus_1);
  73.  
  74. %
  75. % add code to model periodic red light
  76. %
  77. %if ((mod(ts(n), 2/60) < 1/60) && i == Nx0)
  78. % if i == 200
  79. % if ts(n) < 1/60 || (2/60 <= ts(n) && ts(n) < 3/60)
  80. % flux_plus = 0;
  81. % end
  82. % end
  83.  
  84. %
  85. % compute the density
  86. %
  87. rho(n,i)=rho(n-1,i)+delt_by_delx*(flux_plus-flux_minus);
  88.  
  89. % save flux for the previous x
  90. flux_minus = flux_plus;
  91.  
  92. end
  93.  
  94. %
  95. % compute the sample trajectories, fix the code below!
  96. %
  97. for j=1:Ntraj
  98. traj(n,j) = traj(n-1,j) + delta_t*u_max*(1-interpolate_rho(traj(n-1,j), rho(n-1,:), xs) / rho_max);
  99. end
  100.  
  101. end
  102.  
  103. [XX YY] = meshgrid(xs,ts);
  104. pcolor(XX,YY,rho);
  105. shading('interp');
  106.  
  107. hold on;
  108. plot(traj(:,1), ts', 'w-')
  109. plot(traj(:,2), ts', 'w-')
  110. plot(traj(:,3), ts', 'w-')
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement