生成有限差分的矩阵

Question

我正在为 2D PDE 问题实施有限差分方案。我希望避免使用循环来生成有限差分。例如，要生成 u(x,y)_xx 的二阶中心差分，我可以将 u(x,y) 乘以以下内容：

是否有一个很好的矩阵表示 u_xy = (u_{i+1,j+1} + u_{i-1,j-1} - u_{i-1,j+1 } - u_{i+1,j-1})/(4dxdy)?这是一个更难编码的问题，因为它是二维的——我想用 u(x,y) 乘以一些矩阵以避免循环。非常感谢！

Answer 1

显然您可以自己创建矩阵，但在 Matlab 中有 tridiag 用于此目的。

例如

>> full(gallery('tridiag',5,-1,2,-1))

答案=

 2    -1     0     0     0
-1     2    -1     0     0
 0    -1     2    -1     0
 0     0    -1     2    -1
 0     0     0    -1     2

Answer 2

如果你的点存储在 N-by-N 矩阵中，那么，正如你所说，左乘你的有限差分矩阵给出关于 u_{xx} 的二阶导数的近似值。右乘有限差分矩阵的转置相当于一个近似u_{yy}。您可以通过左乘和右乘得到混合导数 u_{xy} 的近似值，例如中心差分矩阵

delta_2x =

     0     1     0     0     0
    -1     0     1     0     0
     0    -1     0     1     0
     0     0    -1     0     1
     0     0     0    -1     0

（然后除以系数 4*Dx*Dy），所以类似

U_xy = 1/(4*Dx*Dy) * delta_2x * U_matrix * delta_2x';

如果将矩阵转换为 N^2 向量

U_vec = U_matrix(:);

那么这些运算符可以用Kronecker product, implemented in MATLAB as kron表示：我们有

A*X*B = kron(B',A)*X(:);

所以对于你的有限差分矩阵

U_xy_vec = 1/(4*Dx*Dy)*(kron(delta_2x,delta_2x)*U_vec);

如果您有一个 N-by-M 矩阵 U_mat，则左矩阵乘法等于 kron(eye(M),delta_2x_N)，右乘等于 kron(delta_2y_M,eye(N))，其中 delta_2y_M (delta_2x_N)是M-by-M(N-by-N)中心差分矩阵，所以运算是

U_xy_vec = 1/(4*Dx*Dy) * kron(delta_2y_M,delta_2y_N)*U_vec;

这是一个 MATLAB 代码示例：

N = 20;
M = 30;
Dx = 1/N;
Dy = 1/M;

[Y,X] = meshgrid((1:(M))./(M+1),(1:(N))/(N+1));

% Example solution and mixed derivative (chosen for 0 BCs)
U_mat = sin(2*pi*X).*(sin(2*pi*Y.^2));
U_xy = 8*pi^2*Y.*cos(2*pi*X).*cos(2*pi*Y.^2);

% Centred finite difference matrices
delta_x_N = 1/(2*Dx)*(diag(ones(N-1,1),1) - diag(ones(N-1,1),-1));
delta_y_M = 1/(2*Dy)*(diag(ones(M-1,1),1) - diag(ones(M-1,1),-1));

% Cast U as a vector
U_vec = U_mat(:);

% Mixed derivative operator
A = kron(delta_y_M,delta_x_N);

U_xy_num = A*U_vec;
U_xy_matrix = reshape(U_xy_num,N,M);

subplot(1,2,1)
contourf(X,Y,U_xy_matrix)
colorbar
title 'Numeric U_{xy}'
subplot(1,2,2)
contourf(X,Y,U_xy)
colorbar
title 'Analytic U_{xy}'

Answer 3

使用 MATLAB 中可用的稀疏函数生成有限差分近似矩阵是一个不错的选择。它节省了大量（实际上非常多）内存...

生成有限差分的矩阵

Matrix to generate finite difference

matlab

numerical-methods

pde

numerical-analysis