Диффузия

Три разных подхода к моделированию диффузии. Рассматриваем одномерное уравнение.

Метод конечных разностей

% Решает одномерное уравнение диффузии на прЯмоугольнике,
% с начальным условием u(t(1),x)= f
% и граничными условиЯми u(t,х(1)) = ua, u(t,х(end)) = ub
D = 2;
xmin = -5; xmax = 5;
tmin = 0; tmax = 4;
J = 10; % число узлов по x
N = 20; % число узлов по t
f = zeros(1,J); f( round(J/2) ) = 50;
ua = 0; ub = 0;
dx = (xmax-xmin)/J;
dt = (tmax-tmin)/N;
s = D*dt/dx^2;
u = zeros(N,J);
u(1,:) = f;
for n = 1:N-1
    u(n+1,2:J-1) = s*(u(n,3:J) + u(n,1:J-2)) + (1-2*s)*u(n,2:J-1);
    u(n+1,1) = ua;
    u(n+1,J) = ub;
end
surf(u)

Однако не все так просто! Численное решение дифференциальных уравнений в частных производных чревато «опасностями», например, если мы зададим N = 10, то получим явно неправдоподобный результат, т.е столкнемся с т.н. неустойчивостью численного метода. Но подробности оставим курсу «Численные методы»...

Метод клеточных автоматов

А здесь даже уравнение писать не нужно. В первый момент времени примесь есть только в центральной ячейке. Концентрация примеси в ячейке на следующем шаге частью остается в самой ячейке, частью переходит в соседние — правую и левую. Если течения нет, то эти доли одинаковы - по 1/3.

ncells = 25; % количество Ячеек
N = 20; % количество шагов по времени

L = 1/3; % влево
S = 1/3; % на месте
R = 1/3; % вправо

c = zeros(ncells,1); % Массив значений концентрации
cmax = 100;          % Максимум концентрации
c(ncells/2) = cmax;   % Начальное значение концентрации
x = 2:ncells-1;

for i=1:N   % цикл по времени
    plot(c) % рисуем текущие значениЯ концентрации          
    c(x) = c(x-1)*L + c(x)*S + c(x+1)*R; % считаем новые
    axis([1 ncells 0 cmax]), grid on % диапазоны по осЯм
    xlabel('Position')      
    ylabel('Value')
    title(['After Time Step: ' num2str(i)])
    drawnow   % Рисует график немедленно
                     % Обычно, MATLAB строит график только в конце программы
end

Мы можем усовершенствовать нашу модель, если введем течение. Тогда доли L, S и R изменятся. Например, если течение направлено вправо, то они могут быть 1/10, 3/10 и 6/10 соответственно. В дифференциальном уравнении этого можно было бы добиться, делая переменным коэффициент диффузии. Правда тогда придется переделать и разностное уравнение.

Заметьте, что концы трубы в программе исключены из рассмотрения (x = 2:ncells-1;). Чтобы рассмотреть их, нужно с помощью условного оператора проверять, достигнут ли конец трубы. Тогда можно сделать концы, например, отражающими, т.е. частицы будут отскакивать от них вглубь трубы.

Источник можно задавать разовой концентрации (как в примере), а можно постоянной. Можно задавать перемещение источника вдоль трубы.

Модель очевидно обобщается на 2-х и 3-х мерный случаи, только количество соседей станет больше. Тогда в ней можно вводить препятствия различных видов (модель распространения загрязнения с домами и лесом).

Вспомнив, про то, что движение каждой частицы примеси происходит случайным образом, мы можем применить для моделирования ее перемещений...

Метод Монте-Карло

Частица может на каждом шаге двинуться вправо, влево или остаться на месте. Через некоторое время (количество шагов) она окажется на каком-то расстоянии от исходной точки. Если запустить множество частиц и позволить им всем двигаться по одинаковым законам (вправо, влево или стоять), а через некоторое время посмотреть, сколько из них удалились на то или иное расстояние от источника, то мы и получим распределение концентраций (количество части в единице объема).

n = 20; % число Ячеек
m = 500; % число частиц
x = (n/2)*ones(1,m);
dice = [-1,0,+1];

for k = 1:20
    dx = dice( round(1+2*rand(1,m)) );
    hist(x,m)
    axis([0,n,0,m/2]); axis square; grid on
    x = x + dx;
    drawnow
end