对于单值多元函数，有些函数难以获得解析导数，因此编程实现了导数的数值算法，基本思想就是用一阶差分代替导数：

其中${\mathbf{x}}_0$表示初值，$d\mathbf{x}$表示$\mathbf{x}$在某一方向上的微小增量，设增量大小$dx=10^{-6}$，根据上式得出该方向上的导数，求解所有方向上的导数，形成数量场对向量的一阶导数（梯度）$\mathbf{g}$。

function [g]=grad_numerical(f,x0,dx)
    %计算数量场f在x0处的梯度g（数值计算，差分距离dx）仅适用于实数
    n=length(x0);
    g=zeros(n,1);
    for i=1:n
        x1=x0;
        x1(i)=x1(i)+dx;
        g(i)=(f(x1)-f(x0))/dx;%偏导定义
    end
end

如果函数是实值函数，而自变量是复数，那么可以分别对实部和虚部求导，相应代码修改为：

function [g]=grad_numerical_CtoR(f,x0,dx)
    %计算数量场f在x0处的梯度g（数值计算，差分距离dx）实值函数对复数求导
    n=length(x0);
    %% 实数的梯度
    gr=zeros(n,1);
    for i=1:n
        x1=x0;
        x1(i)=x1(i)+dx;
        gr(i)=(f(x1)-f(x0))/dx;%偏导定义
    end
    %% 复数的梯度
    gi=zeros(n,1);
    for i=1:n
        x1=x0;
        x1(i)=x1(i)+1i*dx;
        gi(i)=(f(x1)-f(x0))/dx;%偏导定义
    end
    g=gr+1i.*gi;%合成梯度
end

（按理说实值函数不满足柯西黎曼条件，不存在导数，也不知道我这么做对不对）

先把实数的浅浅验证一下：

fx =@(x) norm(x,2)^2; %目标函数（标量）
gx=@(x) 2.*x;%解析梯度
rng(0)
x0=rand(3,1);
gx_num=grad_numerical(fx,x0,1e-6)%数值计算梯度
gx_ana=gx(x0)%梯度解析式

输出结果一样（近似）

专栏里还有数值二阶导数（海森矩阵）的计算方法。