我尝试使用有限元求解二维扩散方程:
numx = 101; % number of grid points in x
numy = 101;
numt = 1001; % number of time steps to be iterated over
dx = 1/(numx - 1);
dy = 1/(numy - 1);
dt = 1/(10*(numt - 1));
x = 0:dx:1; %vector of x values, to be used for plotting
y = 0:dy:1; %vector of y values, to be used for plotting
P = zeros(numx,numy,numt); %initialize everything to zero
% Initial Condition
mu = 0.5;
sigma = 0.05;
for i=2:numx-1
for j=2:numy-1
P(i,j,1) = exp(-( (x(i)-mu)^2/(2*sigma^2) + (y(j)-mu)^2/(2*sigma^2))/(2*sigma^2)) / sqrt(2*pi*sigma^2);
end
end
% Diffusion Equation
for k=1:numt
for i=2:numx-1
for j=2:numy-1
P(i,j,k+1) = P(i,j,k) + .5*(dt/dx^2)*( P(i+1,j,k) - 2*P(i,j,k) + P(i-1,j,k) ) + .5*(dt/dy^2)*( P(i,j+1,k) - 2*P(i,j,k) + P(i,j-1,k) );
end
end
end
figure(1)
surf(P(:,:,30))
P的初始条件是高斯分布。
但是,当增加k(在我的代码中k = 30)时,P会爆炸。但是,P应该减小,因为它是扩散方程的解。
我不知道P(i,j,k + 1)中的问题,如何解决此问题?
谢谢
P(1, :, :)
,P(101, :, :)
,P(:, 1, :)
和P(:, 101 , :)
的值在您的代码中始终为零。这意味着在二维框的边界处P始终为零。它不遵守扩散定律。
首先,应该给边界定一个初始条件。
% Initial Condition
mu = 0.5;
sigma = 0.05;
for i=1:numx
for j=1:numy
P(i,j,1) = exp(-( (x(i)-mu)^2/(2*sigma^2) + (y(j)-mu)^2/(2*sigma^2))/(2*sigma^2)) / sqrt(2*pi*sigma^2);
end
end
然后在扩散过程中计算边界处的P值。
% Diffusion Equation
for k=1:numt
for i=2:numx-1
for j=2:numy-1
P(i,j,k+1) = P(i,j,k) + .5*(dt/dx^2)*( P(i+1,j,k) - 2*P(i,j,k) + P(i-1,j,k) ) + .5*(dt/dy^2)*( P(i,j+1,k) - 2*P(i,j,k) + P(i,j-1,k) );
end
end
% calculate the value of P at the boudary
P(1, :, k + 1) = 2 * P(2, :, k + 1) - P(3, : , K + 1);
P(numx, :, k + 1) = 2 * P(numx - 1, :, k + 1) - P(numx - 2, : , K + 1);
P(:, 1, k + 1) = 2 * P(:, 2, k + 1) - P(:, 3 , K + 1);
P(:, numy, k + 1) = 2 * P(:, numy - 1, k + 1) - P(:, numy - 2 , K + 1);
end