求切平面
我们将展示如何使用有限差分近似计算函数的梯度，并在给定点处构造一个切平面到曲面。
In [ ]:
Pkg.add( "Zygote" )
让我们设置初始函数
让我们创建一个表面函数进行进一步分析:
In [ ]:
f( xy ) = xy[1].^2 .+ xy[2].^2;
在这个例子中，对我们来说，所有参数都使用向量传递给函数是很重要的
要获得感兴趣点处切平面的参数，我们需要对此函数取导数。
在Julia中，有几种方法可以获得函数的导数，例如，使用自动微分和库 Zygote.
In [ ]:
using Zygote
fx, fy = f'( [-1, 3] )
该库可以区分Julia语言中非常广泛的函数：例如，它们的可微分代码可以包含复数，循环，递归，条件操作和来自标准库（例如，FFTW）的命令调用。
点处切平面的方程 你可以设置如下:
系数 和 我们通过使用 f'( 1,2 ). 设置功能 ∇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]:
数值解(diff函数)
现在让我们假设我们只有数据，而不是函数，我们需要两个坐标中的偏导数。 我们得到某个坐标网格的每个点上的导数的近似值:
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 )
Out[0]:
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=(-88, 8.5) )
Out[0]:
结论
我们以两种不同的方式构建了函数的切平面:
-通过用Julia语言区分函数代码（在任何时候）,
-使用数值微分（在某个坐标网格上）。