寻找切平面¶
我们将演示如何利用有限差分近似计算函数的梯度，并在给定点处构建曲面的切平面。
In [ ]:
Pkg.add(["Zygote"])
In [ ]:
Pkg.add( "Zygote" )
让我们定义原始函数¶
让我们创建一个曲面函数，以便进一步分析：
$$f(x,y) = x^2 + y^2$$
In [ ]:
f( xy ) = xy[1].^2 .+ xy[2].^2;
在本例中，重要的是所有参数都要通过向量
要获得切平面在兴趣点处的参数，我们需要对该函数进行导数运算。
在 Julia 中，有多种获取函数导数的方法，例如使用自求导数和库Zygote 。
In [ ]:
using Zygote
fx, fy = f'( [-1, 3] )
Out[0]:
这个库可以微分 Julia 中非常广泛的函数：例如，其微分代码可以包含复数、循环、递归、条件运算以及调用标准库（如 FFTW）的命令。
切平面在点$P=[x_0,y_0]$ 处的方程可以如下指定：
$$z = f(x_0, y_0) + \frac{\partial f(x_0,y_0)}{\partial x} (x - x_0) + \frac{\partial f(x_0,y_0)}{\partial y} (y - y_0)$$
我们利用f'( 1,2 ) 得到了系数$\frac{\partial x}{\partial y}$ 和$\frac{\partial x}{\partial y}$ 。让我们定义函数∇f(x,y) ，它可以让我们找到任意点的切平面值：
In [ ]:
x₀,y₀ = -1,3
fx,fy = f'( [x₀,y₀] )
∇f(x,y) = f( [x₀,y₀] ) .+ fx.*(x.-x₀) .+ fy.*(y.-y₀)
Out[0]:
让我们推导出原始函数f(x,y) 和在点P = [1,2] 上的切平面：
In [ ]:
gr()
x_range = -5:0.25:5
y_range = -5:0.25:5
xx,yy = repeat( x_range, inner = length(y_range)), repeat( y_range, outer = length(x_range));
Nx = Ny = round( Int32, sqrt(length(xx)) );
xm = reshape( xx, Nx, Nx );
ym = reshape( yy, Ny, Ny );
wireframe( x_range, y_range, ∇f(xm,ym), cbar=false )
surface!( x_range, y_range, f([xm, ym]), c=:viridis, cbar=false )
scatter!( [x₀], [y₀], [f([x₀, y₀])], c=:white, leg=false )
Out[0]:
让我们从另一个角度来看这个图形。
In [ ]:
plot!( camera=(-86, 9) )
Out[0]:
数值解法（差分函数）¶
现在假设我们只有数据，没有函数，需要两个坐标上的部分导数。让我们求出某个坐标网格上每一点的导数近似值：
In [ ]:
dfx = diff( f( [xm, ym] ), dims=2 );
dfy = diff( f( [xm, ym] ), dims=1 );
plot( surface(dfx), surface(dfy), cbar=false, size=(500,250) )
Out[0]:
要在(x₀, y₀) （坐标网格上的其中一点）处绘制切线，我们需要找到该点的索引：
In [ ]:
t = (xm .== x₀) .& (ym .== y₀)
indt = findfirst( t )
In [ ]:
dT = yy[2]-yy[1] # Расстояние между всеми точками: 0.25
fx0 = dfx[ indt ] ./ dT;
fy0 = dfy[ indt ] ./ dT;
同样，我们会得到相同的切线图形：
In [ ]:
∇f2( x,y ) = f(x₀,y₀) .+ fx0.*(x.-x₀) .+ fy0.*(y.-y₀)
wireframe( x_range, y_range, ∇f2(xm,ym) )
surface!( x_range, y_range, f([xm, ym]), c=:viridis, cbar=false )
plot!( camera=(-86, 9) )
Out[0]:
结论¶
我们用两种不同的方法构造了函数的切平面：
- 在 Julia 中对函数代码*进行微分（在任意点）、
- 数值微分（在某个坐标网格上）。