0

我正在求解泊松方程，并想绘制精确解的误差与网格点数的关系。我的代码是：

function [Ntot,err] = poisson(N)


nx = N;                % Number of steps in space(x)
ny = N;                % Number of steps in space(y)       
Ntot = nx*ny;
niter = 1000;           % Number of iterations 
dx = 2/(nx-1);          % Width of space step(x)
dy = 2/(ny-1);          % Width of space step(y)
x = -1:dx:1;             % Range of x(-1,1) 
y = -1:dy:1;             % Range of y(-1,1) 
b = zeros(nx,ny);       
dn = zeros(nx,ny);      


% Initial Conditions
d = zeros(nx,ny);                  
u = zeros(nx,ny);

% Boundary conditions
d(:,1) = 0;
d(:,ny) = 0;
d(1,:) = 0;                  
d(nx,:) = 0;


% Source term
b(round(ny/4),round(nx/4)) = 3000;
b(round(ny*3/4),round(nx*3/4)) = -3000;


i = 2:nx-1;
j = 2:ny-1;

% 5-point difference (Explicit)
for it = 1:niter
    dn = d;
    d(i,j) = ((dy^2*(dn(i + 1,j) + dn(i - 1,j))) + (dx^2*(dn(i,j + 1) + dn(i,j - 1))) - (b(i,j)*dx^2*dy*2))/(2*(dx^2 + dy^2));
    u(i,j) = 2*pi*pi*sin(pi*i).*sin(pi*j);

    % Boundary conditions 
    d(:,1) = 0;
    d(:,ny) = 0;
    d(1,:) = 0;                  
    d(nx,:) = 0;
end


%     
% 
% err = abs(u - d);

我得到的错误是：

下标分配尺寸不匹配。

泊松误差（第 39 行）

u(i,j) = 2*pi*pi*sin(pi*i).*sin(pi*j);

我不确定为什么它不在每个网格点计算 u。我尝试将其从 for 循环中取出，但这并没有帮助。任何想法，将不胜感激。

4

1 回答 1

0

这是因为i和j都是 1×(N-2) 向量，所以u(i, j)也是一个 (N-2)×(N-2)矩阵。但是，表达式2*pi*pi*sin(pi*i).*sin(pi*j)是一个 1×(N-2)向量。
尺寸显然不匹配，因此出现错误。

我不确定，但我猜您打算执行以下操作：

u(i,j) = 2 * pi * pi * bsxfun(@times, sin(pi * i), sin(pi * j)');

或者，您可以使用基本矩阵乘法来生成 (N-2)-by-(N-2)，如下所示：

u(i, j) = 2 * pi * pi * sin(pi * i') * sin(pi * j); %// Note the transpose

PS：建议不要使用“i”和“j”作为变量名。

于 2013-03-11T23:30:10.600 回答