我正在尝试通过 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
我的共轭梯度求解器可以单独使用,但我之前没有将它用作函数。任何帮助将不胜感激。