如何解决如何向量化这个五点差分代码以快速计算矩阵的导数?
我正在尝试使用 five point method 计算二维矩阵中每个维度的偏一阶导数,即 dF/dx 和 dF/dy。我通过遍历点成功地做到了这一点:
dF_dx = zeros(size(F));
dF_dy = zeros(size(F));
% Derivative with respect to y for each x value (Apply to all columns simultaneously)
dF_dy(2,1:(Nx-1)) = ( F(3,1:(Nx-1)) - F(1,1:(Nx-1)) )/(2*dy);
for m = 3:(Ny-2)
dF_dy(m,1:(Nx-1)) = ( F(m-2,1:(Nx-1)) - F(m+2,1:(Nx-1)) + 8*F(m+1,1:(Nx-1)) - 8*F(m-1,1:(Nx-1)) )/(12*dy);
end
dF_dy(Ny-1,1:(Nx-1)) = ( F(Ny,1:(Nx-1)) - F(Ny-2,1:(Nx-1)) )/(2*dy);
% Derivative with respect to x for each y value (Apply to all rows simultaneously)
dF_dx(2:(Ny-1),2) = ( F(2:(Ny-1),3) - F(2:(Ny-1),1) )/(2*dx);
for n = 3:(Nx-2)
dF_dx(2:(Ny-1),n) = ( F(2:(Ny-1),n-2) - F(2:(Ny-1),n+2) + 8*F(2:(Ny-1),n+1) - 8*F(2:(Ny-1),n-1) )/(12*dx);
end
dF_dx(2:(Ny-1),(Nx-1)) = ( F(2:(Ny-1),Nx) - F(2:(Ny-1),Nx-2) )/(2*dx);
是否有一种聪明的 Matlab 方法可以将其矢量化以使其执行速度更快?我已经阅读了几个似乎可以解决问题的函数(例如 diff()、circshift() 甚至 kron() ),但我不确定如何使用它们来解决这个问题问题。
(我的计划是稍后将其实现为 gpuArray,如果这与解决方案相关)。
谢谢!
第一次编辑
通过查看 gradient() 的源代码,我能够制作以下矢量化版本(即没有循环):
F = rand(5,8);
Nx = size(F,2);
Ny = size(F,1);
dx = 2;
dy = 3;
dF_dx = zeros(size(F));
dF_dy = zeros(size(F));
dF_dx(2:(Ny-1),3:Nx-2) = (F(2:(Ny-1),1:Nx-4) - F(2:(Ny-1),5:Nx) + 8*F(2:(Ny-1),4:Nx-1) - 8*F(2:(Ny-1),2:Nx-3))/(12*dx);
dF_dx(2:(Ny-1),1) )/(2*dx);
dF_dx(2:(Ny-1),Nx-2) )/(2*dx);
dF_dy(3:Ny-2,1:(Nx-1)) = (F(1:Ny-4,1:(Nx-1)) - F(5:Ny,1:(Nx-1)) + 8*F(4:Ny-1,1:(Nx-1)) - 8*F(2:Ny-3,1:(Nx-1)))/(12*dy);
dF_dy(2,1:(Nx-1)) )/(2*dy);
dF_dy((Ny-1),1:(Nx-1)) )/(2*dy);
有没有办法在没有所有索引的情况下做到这一点(我怀疑这是最耗时的)?
第二次编辑
我现在已经实施了 Cris 和 chtz 的建议,并使用 conv2()
的卷积实现了这一点,如下所示:
F = rand(500,600);
Nx = size(F,1);
dx = 2;
dy = 3;
[dF_dx,dF_dy] = partial_derivatives(F,dx,dy,Nx,Ny)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dF_dx,Ny)
kernx = [-1 8 0 -8 1]/(12*dx); % Convolution kernel for x dimension
kerny = [-1;8;0;-8;1]/(12*dy); % Convolution kernel for y dimension
%%%%%%%% Partial derivative across x dimension %%%%%%%%
dF_dx = conv2(F,kernx,'same') ; % Internal mesh points (five-point method)
% Second and penultimate mesh points (two-point method)
dF_dx(2:(Ny-1),Nx-2) )/(2*dx);
dF_dx(:,[1 Nx]) = 0; % Set boundary conditions
dF_dx([1 Ny],:) = 0; %
%%%%%%%% Partial derivative across x dimension %%%%%%%%
dF_dy = conv2(F,kerny,'same') ; % Internal mesh points (five-point method)
% Second and penultimate mesh points (two-point method)
dF_dy(2,1:(Nx-1)) )/(2*dy);
dF_dy(:,Nx) = 0; % Set boundary conditions
dF_dy([1 Ny],:) = 0; %
end
作为一个 gpuArray,与第 1 次编辑中的矢量化版本相比,这并没有提供任何显着的改进。有没有明显的改进方法?谢谢。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。