Gradient.m

% 计算三维标量场 U(x, y, z) 的梯度分量

function [DxU, DyU, DzU] = Gradient(x, y, z, U)
[Nx, Ny, Nz] = size(U);

% 初始化输出矩阵
DxU = zeros(Nx, Ny, Nz);
DyU = zeros(Nx, Ny, Nz);
DzU = zeros(Nx, Ny, Nz);

% 计算x方向的梯度
DxU(2:Nx-1,:,:) = (U(3:Nx,:,:) - U(1:Nx-2,:,:))./(x(3:Nx) - x(1:Nx-2));
DxU(1,:,:) = (-U(3,:,:) + 4*U(2,:,:) - 3*U(1,:,:)) / (-x(3) + 4*x(2) - 3*x(1));
DxU(Nx,:,:) = (3*U(Nx,:,:) - 4*U(Nx-1,:,:) + U(Nx-2,:,:)) / (3*x(Nx) - 4*x(Nx-1) + x(Nx-2));

% 计算y方向的梯度
DyU(:,2:Ny-1,:) = (U(:,3:Ny,:) - U(:,1:Ny-2,:))./(y(3:Ny) - y(1:Ny-2))';
DyU(:,1,:) = (-U(:,3,:) + 4*U(:,2,:) - 3*U(:,1,:)) / (-y(3) + 4*y(2) - 3*y(1));
DyU(:,Ny,:) = (3*U(:,Ny,:) - 4*U(:,Ny-1,:) + U(:,Ny-2,:)) / (3*y(Ny) - 4*y(Ny-1) + y(Ny-2));

% 计算z方向的梯度
DzU(:,:,2:Nz-1) = permute( ...
    permute(U(:,:,3:Nz) - U(:,:,1:Nz-2), [3 2 1]) ...
    ./(z(3:Nz) - z(1:Nz-2)), [3 2 1]);
DzU(:,:,1) = (-U(:,:,3) + 4*U(:,:,2) - 3*U(:,:,1)) / (-z(3) + 4*z(2) - 3*z(1));
DzU(:,:,Nz) = (3*U(:,:,Nz) - 4*U(:,:,Nz-1) + U(:,:,Nz-2)) / (3*z(Nz) - 4*z(Nz-1) + z(Nz-2));
end