Finding the tangent plane

We will show how to approximately calculate the gradient of a function using finite differences and construct a tangent plane to the surface at a given point.

Pkg.add(["Zygote"])

Pkg.add( "Zygote" )

Let us define the original function

Let's create a surface function for further analysis:

$$f(x,y) = x^2 + y^2$$

f( xy ) = xy[1].^2 .+ xy[2].^2;

In this example, it is important that all parameters are passed to the function using the vector

To obtain the parameters of the tangent plane at the point of interest, we need to take the derivative of this function.

In Julia, there are several ways to get the derivative of functions, for example, using autodifferentiation and the library Zygote.

using Zygote
fx, fy = f'( [-1, 3] )

2-element Vector{Float64}:
 -2.0
  6.0

This library can differentiate a very wide range of functions in Julia: for example, their differentiable code can contain complex numbers, loops, recursion, conditional operations, and calls to commands from standard libraries (e.g. FFTW).

The equation of the tangent plane at the point $P=[x_0,y_0]$ can be specified as follows:

$$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)$$

We obtained the coefficients $\frac{\partial x}{\partial y}$ and $\frac{\partial x}{\partial y}$ using f'( 1,2 ). Let us define the function ∇f(x,y), which allows us to find the values of the tangent plane at any point:

x₀,y₀ = -1,3
fx,fy = f'( [x₀,y₀] )
∇f(x,y) = f( [x₀,y₀] ) .+ fx.*(x.-x₀) .+  fy.*(y.-y₀)

∇f (generic function with 1 method)

Let's derive the original function f(x,y) and the tangent plane at the point P = [1,2]:

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 )

Let's look at this graph from another angle.

plot!( camera=(-86, 9) )

Numerical solution (diff function)

Now suppose we only have the data, not the function, and we need partial derivatives at both coordinates. Let us get an approximate value of the derivatives at each point of some coordinate grid:

dfx = diff( f( [xm, ym] ), dims=2 );
dfy = diff( f( [xm, ym] ), dims=1 );

plot( surface(dfx), surface(dfy), cbar=false, size=(500,250) )

To plot the tangent at (x₀, y₀), one of the points on the coordinate grid, we find the index of that point:

t = (xm .== x₀) .& (ym .== y₀)
indt = findfirst( t )

dT = yy[2]-yy[1] # Расстояние между всеми точками: 0.25
fx0 = dfx[ indt ] ./ dT;
fy0 = dfy[ indt ] ./ dT;

Again we get the same graph of the tangent:

∇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) )

Conclusion

We have constructed the tangent plane to a function in two different ways:

• by differentiation of the function code in Julia (at any point),
• by numerical differentiation (on some coordinate grid).