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）x（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 に答える