Двумерная интерполяция вне проекции

Question

Допустим, у меня есть 3d-объект на сетке V = V(a, b, c). Я хочу интерполировать V(a, b + alpha*d, c+d).

Другими словами, определите f(d) = V(a, b + alpha*d, c+d). Я хочу приблизить f. Важно отметить, что затем я хочу применить optimize.root к приближению, поэтому я ценю эффективные вычисления f.

Например,

gamma = 0.5
aGrid = np.linspace(5, 10, 30)
bGrid = np.linspace(4, 7, 40)
cGrid = np.linspace(0.1, 0.5, 20)
A, B, C = np.meshgrid(aGrid, bGrid, cGrid, indexing='ij')
V = A**2 + B*C
# define initial a, b, c
idx = (7, 8, 9)
a, b, c = A[idx], B[idx], C[idx]
# so V(a, b, c) = V[idx]

Наивный подход будет

g = scipy.interpolate.interp2d(bGrid, cGrid, V[7, ...])
f = lambda x: g(b + gamma*x, c + x)

и моя конечная цель:

constant = 10
err = lambda x: f(x) - constant
scipy.optimize.root(err, np.array([5]))

Однако все это выглядит очень грязно и неэффективно. Есть ли более питонный способ сделать это?

Answer 1

Я изменил обозначения, чтобы помочь мне в понимании вопроса (я привык к обозначениям физики).В трехмерном пространстве есть скалярное поле V(x, y, z).Мы определяем параметрическую линию в этом трехмерном пространстве:

f_{x0, y0, z0, v_x, v_y, v_z}(t) = (x0 + v_x*t, y0 + v_y*t, z0 + v_z*t)

Это можно рассматривать как траекторию частицы, начинающейся в точке (x0, y0, z0) и движущейся по прямой линии с вектором скорости (v_x, v_y, v_z).

Мы ищем время t1, чтобы V( f(t1) ) было равно заданному значению V0.Это заданный вопрос?

import numpy as np
from scipy.interpolate import RegularGridInterpolator
from scipy.optimize import root_scalar
import matplotlib.pylab as plt

# build the field
aGrid = np.linspace(5, 10, 30)
bGrid = np.linspace(4, 7, 40)
cGrid = np.linspace(0.1, 0.5, 20)
A, B, C = np.meshgrid(aGrid, bGrid, cGrid, indexing='ij')
V = A**2 + B*C

# Build a continuous field by linear interpolation of the gridded data:
V_interpolated = RegularGridInterpolator((aGrid, bGrid, cGrid), V,
                                         bounds_error=False, fill_value=None)

# define the parametric line
idx = (7, 8, 9)
x0, y0, z0 = A[idx], B[idx], C[idx]
alpha = 0.5
v_x, v_y, v_z = 0, alpha, 1 

def line(t):
    xyz = (x0 + v_x*t, y0 + v_y*t, z0 + v_z*t) 
    return xyz

# Plot V(x,y,z) along this line (to check, is there a unique solution?)
t_span = np.linspace(0, 10, 23)
V_along_the_line = V_interpolated( line(t_span) )

plt.plot(t_span, V_along_the_line);
plt.xlabel('t'); plt.ylabel('V');

# Find t such that V( f(t) ) == V1
V1 = 80 
sol = root_scalar(lambda s: V_interpolated( line(s) ) - V1,
                  x0=.0, x1=10)

print(sol)

#      converged: True
#           flag: 'converged'
# function_calls: 8
#     iterations: 7
#           root: 5.385594973846983

# Get the coordinates of the solution point:
print("(x,y,z)_sol = ", line(sol.root))

# (6.206896551724138, 7.308182102308106, 5.675068658057509)