矩陣生成有限差異

我正在實施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

矩陣生成有限差異

回答

相關問題