对于单值多元函数,有些函数难以获得解析导数,因此编程实现了导数的数值算法,基本思想就是用一阶差分代替导数:
其中
x
0
\mathbf{x}_0
x0表示初值,
d
x
d\mathbf{x}
dx表示
x
\mathbf{x}
x在某一方向上的微小增量,设增量大小
d
x
=
1
0
−
6
dx=10^{-6}
dx=10−6,根据上式得出该方向上的导数,求解所有方向上的导数,形成数量场对向量的一阶导数(梯度)
g
\mathbf{g}
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)%梯度解析式
输出结果一样(近似)
专栏里还有数值二阶导数(海森矩阵)的计算方法。