计算科学 - 继续解决− ∇ ⋅ ( k ( x , y) ∇ u ) = f−∇⋅(k(x,y)∇u)=f - 吾爱随笔录

继续解决− ∇ ⋅ ( k ( x , y) ∇ u ) = f−∇⋅(k(x,y)∇u)=f

计算科学 pde matlab 有限差分

2021-12-26 01:49:28

我正在尝试解决以下问题，我之前已经打开了另一个关于实现的讨论，而且看起来效果很好，可以在这里找到。

对于 n = 500,1000,5000 和 10,000 点的网格，我需要计算具有不完全 Cholesky 和 Jacobi 的预条件共轭梯度解。但对于个人电脑来说，这些是很多要点。正如@Abdullah Ali Sivas 评论的那样：

不需要矩阵，可以利用矩阵的作用来解决泊松问题；特别是如果您实施固定迭代方法，如 Jacobi 或 Gauss-Seidel 方法。

当尺寸巨大，每边有很多点并且我的电脑没有这样的能力时，我希望你能帮助我知道如何获得解决这个问题的方法。

使用句柄功能：

function y=afun(x)
nx = 7;
N = nx*nx;
x1 = linspace(0,1,nx);
y1 = x1;
[X,Y] = ndgrid(x1,y1);
dx = x1(2)-x1(1);
k=@(x,y) 1+x.^2+y.^2;
isDirichlet = (X==0) | (X==1) | (Y==0) | (Y==1);
Ad = zeros(N,5);
Ad(:,1) = -k(X(:),Y(:)+dx/2); 
Ad(:,2) = -k(X(:)+dx*0.5,Y(:)); 
Ad(:,4) = -k(X(:)-dx*0.5,Y(:)); 
Ad(:,5) =  -k( X(:), Y(:)-dx/2 );
Ad(:,3) = -sum( Ad(:,[1,2,4,5]), 2 ); 
idx = find(isDirichlet(:)); 
L=Ad(:,1);L(idx)=[];B=Ad(:,2);B(idx)=[];
C=Ad(:,3);C(idx)=[];U=Ad(:,4);
U(idx)=[];UP=Ad(:,5);UP(idx)=[];
n=sqrt(size(C,1));
l=B(1:end-1);s=U(2:end);
for i=n:n:(n-1)*n
l(i)=0;s(i)=0;
end
y=([0;l]+[C]+[s;0]+[diag(zeros(n));L(1:(n-1)*n)]+[UP(n+1:end);diag(zeros(n))]).*x;
end

2个回答

根据您的第一个问题，您有一个稀疏矩阵。这意味着您不需要存储所有 $n\times n$ 你的问题的系数。一种对矩阵向量乘积有用的格式是压缩稀疏行 (CSR)。

在 Jacobi 方法的情况下，您可以重新排列

A x = f,

$Ax = f\, ,$

作为

x^{k + 1} = D^{- 1} (f - (L + U) x),

$x^{k + 1} = D^{-1}(f - (L + U) x)\, ,$

和 $D$ 的对角线 $A$ ，和 $L + U = A - D$ .

同样的处理也适用于 Gauss-Seidel 或共轭梯度方法。

首先，我想问你是否确定当你说 $n=10000$ 你的意思是 $10^4$ 一维点？因为那意味着 $10^8$ 点在域上，解决这样一个简单问题的点太多了。你会在之前达到最小错误 $10^8$ 积分，你不会因为你做的额外工作而获得任何东西。

其次，您实际上并不需要矩阵。你需要它的动作。例如，检查链接https://www.mathworks.com/help/matlab/ref/pcg.html上的“使用函数句柄而不是数字矩阵”示例。主要思想是您创建一个函数，该函数接受一个向量并计算 A*x 而无需构造矩阵 A。在您的情况下，它可能是一个 for 循环。

同样的想法也适用于雅可比。但是，对于 ICHOL，您需要明确地使用矩阵，并且无法绕过它。因此，我的第一点。可能，您预计不会在每个方向上使用超过 500 个点（或总共 250,000 个点）。

其它你可能感兴趣的问题

上一篇SQP 中的过度约束下一篇重叠的一维网格