MATLAB非均匀网格梯度计算方法详解

MATLAB非均匀网格梯度计算方法详解

引用

CSDN

1.

https://m.blog.csdn.net/ambu1230/article/details/138198649

在MATLAB中，gradient函数可以很方便地对均匀网格进行梯度计算，但对于非均匀网格，gradient函数却无法直接求解。本文将介绍如何编写自定义函数来计算三维矩阵的梯度，包括理论推导、代码实现和结果验证。

理论推导

1. 内部网格点

对$a_1$和$a_3$两点分别进行泰勒展开，公式如下：

$$
a_3 = a_2 + \dot{a}_2 \Delta x_2 + \frac{1}{2} \ddot{a}_2 \Delta x_2^2 + O(\Delta x_2^3) \textcircled{1} \
a_1 = a_2 - \dot{a}_2 \Delta x_1 + \frac{1}{2} \ddot{a}_2 \Delta x_1^2 + O(\Delta x_1^3) \textcircled{2}
$$

最终得到：

2. 边界点

代码实现

1D函数

function dydx = calc_grad_1D(x,y)
%% 求解一维数组的梯度
%% input1:一维函数坐标-->x
%% input2：一维函数值-->y
dydx = zeros(1,length(x));
for i = 1:length(x)
    if i>1 && i<length(x)
        deltax1 = x(i)-x(i-1);
        deltax2 = x(i+1)-x(i);
        son = (y(i+1)*deltax1^2-y(i-1)*deltax2^2-y(i)*(deltax1^2-deltax2^2));
        mom = (deltax2*deltax1^2+deltax1*deltax2^2);
        dydx(i) = son/mom;
    elseif i==1
        n = (x(3)-x(1))/(x(2)-x(1));
        son = y(i+2)-y(i+1)*n^2-(1-n^2)*y(i);
        mom = (n-n^2)*(x(i+1)-x(i));
        dydx(i)=son/mom;
    elseif i==length(x)
        n = (x(i)-x(i-2))/(x(i)-x(i-1));
        son = y(i-2)-y(i-1)*n^2-(1-n^2)*y(i);
        mom = (n-n^2)*(x(i)-x(i-1));
        dydx(i)=-son/mom;
    end
end
end

3D矩阵

function [dfdx,dfdy,dfdz] = calc_grad_3D(F,X,Y,Z)
%UNTITLED26 此处提供此函数的摘要
%   此处提供详细说明
nx = size(X,1);ny = size(Y,2);nz = size(Z,3);
dfdx = zeros(nx,ny,nz);dfdy = zeros(nx,ny,nz);dfdz = zeros(nx,ny,nz);
for j = 1:ny
    for k = 1:nz
        dfdx(:,j,k) = calc_grad_1D(X(:,j,k),F(:,j,k));
    end
end
for i = 1:nx
    for k = 1:nz
        dfdy(i,:,k) = calc_grad_1D(Y(i,:,k),F(i,:,k));
    end
end
for i = 1:nx
    for j = 1:ny
        dfdz(i,j,:) = calc_grad_1D(Z(i,j,:),F(i,j,:));
    end
end
end

结果验证

具体案例是求解函数$F = x^2 + y^2 + z^2$在三个方向的梯度：

clc;clear
x = 1:10;y = x;z = x;
[X,Y,Z] = ndgrid(x,y,z);
F = X.^3+Y.^2+Z.^3;
%%
[dFdy,dFdx,dFdz] = gradient(F,Y(1,:,1),X(:,1,1),Z(1,1,:));
%%
[dfdx,dfdy,dfdz] = calc_grad_3D(F,X,Y,Z);
%% 理论解与数值解对比
dfdy_ana = 2.*(Y);
dfdy_ana = reshape(dfdy_ana,1000,1);
dfdy = reshape(dfdy,1000,1);
dFdy = reshape(dFdy,1000,1);
c = abs(dfdy-dfdy_ana);
d = abs(dFdy-dfdy_ana);
plot(c,'-o')
hold on
plot(d,'-o')
%% 绘图设置
axis([0 1000 0 2])
legend('My code','MATLAB gradient')
ylabel('误差')

结果如下：

可以看出，MATLAB里的gradient函数由于在边界上采用一阶差分，因此存在误差，而我们编写的函数在内部点和边界点都采用二阶精度，因此误差为0。