矩阵生成有限差异

我正在实施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）这是一个难以编码的问题，因为它是2D的 - 我想用u（x，y）乘一些矩阵来避免循环。非常感谢！

来源

2016-03-15 Mike Miller

什么是您的编程语言/环境？ – karakfa

@karakfa只是在Matlab –

如果点存储在一个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

，所以像

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

（再由因子4*Dx*Dy分）如果你投了矩阵作为N^2矢量

U_vec = U_matrix(:);

那么这些运营商可以可使用Kronecker product表示，在MATLAB中实现为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_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}'

来源

2016-03-15 15:21:20 Steve

谢谢，这是一个很好的答案。如果u（x1，x2）在矩形域上呢？ –

事情变得棘手一点，但这是可能的。你想在每个方向有不同数量的点是吗？ – Steve

是的，我想概括:)我基本上已经完成了你的建议的第一部分。例如。如果它被分解成6x4网格，我会生成两个tridiagonals，4x4（用于u_xx）和6x6（用于u_yy），它与我在原始文章中使用的结构相同。 u *（4x4）为u_xx，（6x6）* u为u_yy。 –

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

例如

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

ANS =

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

来源

2016-03-15 14:21:57 karakfa

对不起，我应该在我的问题更具体。我知道如何在matlab中创建三对角矩阵，它是我想知道的偏导数矩阵的结构。也许我应该问comp sci？ –

对不起，我没注意，我以为你在描述这个矩阵的结构。我会考虑这个。等式中你的意思是'u_ij = ...'而不是'u_xy'。 – karakfa

由于您所依赖的矩阵索引与行索引或列索引不匹配，我不认为使用简单矩阵乘法是可能的。 – karakfa

使用MATLAB提供稀疏的功能产生有限差分逼近矩阵是一个不错的选择。它节省了大量的内存（确实非常多）。 ..

来源

2017-07-21 06:23:38

矩阵生成有限差异

回答

相关问题