Документация Engee
Notebook

Нахождение касательной плоскости

Открыть пример в Engee

Покажем, как приблизительно рассчитать градиент функции при помощи конечных разностей и построим касательную плоскость к поверхности в заданной точки.

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]:
2-element Vector{Float64}:
 -2.0
  6.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)$$

Коэффициенты $\frac{\partial x}{\partial y}$ и $\frac{\partial x}{\partial y}$ мы получили при помощи 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 (generic function with 1 method)

Выведем исходную функцию 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 )
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 (в любой точке),
  • при помощи численного дифференцирования (на некоторой координатной сетке).