0

我正在尝试通过 Matlab 中的共轭梯度法求解具有分布式数组的泊松方程。我不断收到此错误：

    Index exceeds matrix dimensions.

Error in spdiags (line 103)
         a = [a; i i+d(k) B(i+(m>=n)*d(k),k)];

Error in conjgrad_spmd (line 9)
A = distributed(spdiags(ones(n,1),0,AA));

Error in poisson_conjgrad (line 40)
x = conjgrad_spmd(A, B, n)

我的代码设置矩阵来解决 Ax=b：

  % Solve the Poisson equation via Conjugate Gradient Method

a = -1;                     % Boundary
b = 1;                      % Boundary
n = 50; 
h = (b - a)/(n + 1);        % Step-size
x = linspace(a, b, n + 2);  % # of Grid points x
y = linspace(a, b, n + 2);  % # of Grid points y

[X,Y] = meshgrid(x,y);      % Create 2D Array
X = X';                     % Transpose X, ensure I am getting the interior points (i, j)
Y = Y';                     % Transpose Y, ensure I am getting the interior points (i, j)

Iint = 2:(n + 1);           % Interior x index
Jint = 2:(n + 1);           % Interior y index
Xint = X(Iint, Jint);       % Interior Points
Yint = Y(Iint, Jint);

f = @(x,y) sin(pi*x)*sin(pi*y);         % Our f(x, y)
uexact = (pi^2)*sin(pi*X)*sin(pi*Y);    % Exact solution
u = uexact;                             % Sets full array for u, only need to ensure boundary points are included

% 5-Point Stencil
D2 = f(Xint, Yint);                     % Evaluate f(x, y) @ interior points
D2(:,1) = D2(:,1) - u(Iint, 1)/h^2;     % Ensure that D2 includes boundary terms
D2(:,n) = D2(:,n) - u(Iint, n + 2)/h^2;
D2(1,:) = D2(1,:) - u(1, Jint)/h^2;
D2(n,:) = D2(n,:) - u(n + 2, Jint)/h^2;
B = reshape(D2, n*n, 1);                 % Make D2 into a column vector

% Create A matrix for solving Ax = b
I = speye(n);
e = ones(n, 1);
T = spdiags([e -4*e e],[-1 0 1], n, n);
S = spdiags([e e], [-1 1], n, n);
A = (kron(I,T) + kron(S,I)) / h^2;

% Solve Ax = b for x using Conjugate Gradient Method
tic();
x = conjgrad_spmd(A, B, n)
toc() % Elapsed Time

还有我的共轭梯度求解器：

function [x] = conjgrad_spmd(AA, r, n) % for Ax=b A=AA, b=r

% Conjugate Gradient Experiment. Ax=b
%n = 500; density = 0.5;  tau = 0.01;
%AA = sprandsym(n,density);
AA = -1 + 2*AA;
%flag = abs(AA) > tau;
%AA(flag) = 0;
A = distributed(spdiags(ones(n,1),0,AA));

x = distributed.zeros(n,1);
r = distributed(n,1); % b
spmd
    Niter = 100; 
    err = zeros(Niter,2);
    p = r;
    rsold = dot(r,r);
    for k=1:100
        Ap = A*p; 
        a = rsold / dot(p,Ap);
        x = x + a*p;
        r = r - a*Ap;
        rsnew = dot(r,r);
        sqrtrsnew = sqrt(rsnew);
        err(k,:) = [k sqrtrsnew];
        if sqrtrsnew < eps
            iterstop = k;
            break;
        else
            iterstop = Niter;
        end
        p = r + (rsnew/rsold)*p;
        rsold = rsnew;
    end
end

我的共轭梯度求解器可以单独使用，但我之前没有将它用作函数。任何帮助将不胜感激。

4

1 回答 1

0

错误在 A = Distributed(spdiags(ones(n,1),0,AA)); 尺寸应该是 n*n ： A = Distributed(spdiags(ones(n*n,1),0,AA));

于 2013-05-06T18:18:01.267 回答