How to numerically simulate molecular diffusion

[Ante]

I want to simulate molecular diffusion in MATLAB.

The idea is as follows:

http://www.math.chalmers.se/Math/Grundutb/GU/MVG300/V11/Laborationer/anim1.jpg

The picture represents a gas container and the red part is a barrier which the gas molecules cannot pass through except for the small area above it. The molecules are on the left side at the beginning, and moves RANDOMLY and we will animate this with sequences of pictures in a loop - and after one iteration the molecules position will uppdate.

After 1000 iterations:

http://www.math.chalmers.se/Math/Grundutb/GU/MVG300/V11/Laborationer/anim2.jpg

(x, y) is the initial position of the molecules and the new position is (x + dx, y + dy) where dx = randomnumber, dy = randomnumber.

The particles will hit the barrier and the containers wall, and we will simulate it as follows: If(x + dx, y + dy) is outside the container or within the barrier then we place the particle on the wall. ASay that (x + dx, y + dy) = (-0.21, 0.74).We let the particles position to be (0, 0.74) (analog in other cases).

When the naimation is done, a new window will open and in this the relative amount of particles on the left side will plot as a function of the iterations.

Something like this:

http://www.math.chalmers.se/Math/Grundutb/GU/MVG300/V11/Laborationer/anim3.jpg


Question:
How do I write a MATLAB program diff_anim with above ideas? diff_anim invokes with two parameters diff_anim(no_of_part, no_of_iterations). no_of_part is the amount of particles and no_of_iterations is thje amount of iterations. Particles initial position should randomly generate in the intervall 0 <= x <= 0.9, 0 <= y <= 1.

 

[NateD]

0. % diff_anim.mfunction diff_anim(no_of_part, no_of_iterations) % Initialize the figure figure; % Generate 'no_of_part' random initial positions x = 0.9 * rand(1, no_of_part); y = 1.0 * rand(1, no_of_part); % Iterate the diffusion process for iter = 1:no_of_iterations % Generate random displacements in (x,y) direction dx = 2*rand(1, no_of_part)-1; dy = 2*rand(1, no_of_part)-1; % Move particles x = x + dx; y = y + dy; % Check that particles are not outside the container or in the barrier x(x > 0.9) = 0.9; x(x < 0) = 0; y(y > 1.0) = 1.0; y(y < 0) = 0; % Remove particles in barrier x(y > 0.7) = 0; y(y > 0.7) = 0; % Animate clf; hold on; % Draw container and barrier rectangle('Position', [0, 0, 0.9, 1], 'Facecolor', [0.9 0.9 0.9]); rectangle('Position', [0, 0.7, 0.9, 0.3], 'Facecolor', [1 0 0]); % Plot particles plot(x, y, '.', 'MarkerSize', 10); axis([0 0.9 0 1]); drawnow; % Update particle positions x = x + dx; y = y + dy; end % Compute relative amount of particles on the left side left_side =