Tôi hiện có một khối lượng bao gồm vài triệu hạt phân bố không đều, mỗi hạt có một tính chất (tiềm năng cho những ai tò mò) mà tôi muốn tính lực cục bộ (gia tốc).

np.gradient chỉ hoạt động với dữ liệu cách đều nhau, tôi đang xem nó ở đây:Độ dốc thứ hai trong numpyCần phải nội suy nhưng tôi không thể tìm thấy triển khai spline 3D trong Numpy.

Một số mã sẽ tạo ra dữ liệu đại diện:

nhập numpy dưới dạng np    
từ nhập scipy.spatial cKDTree

x = np.random.uniform(-10, 10, 10000)
y = np.random.uniform(-10, 10, 10000)
z = np.random.uniform(-10, 10, 10000)
phi = np.random.uniform(-10**9, 0, 10000)

kdtree = cKDTree(np.c_[x,y,z])
_, index = kdtree.query([0,0,0], 32) #find 32 hạt gần nhất với gốc tọa độ
#tìm độ dốc tại (0,0,0) bằng cách xem xét 32 hạt gần nhất?  

(Câu hỏi của tôi vàChức năng tính toán độ dốc 3D với các vị trí mẫu cách đều nhauNó tương tự, nhưng dường như không có giải pháp nào nên tôi nghĩ tôi nên hỏi lại. )

Bất kỳ trợ giúp sẽ được đánh giá rất cao.

câu trả lời hay nhất

Đây là cách triển khai Julia thực hiện những gì bạn yêu cầu

sử dụng Hàng xóm gần nhất

n = 3;
k = 32; # để ổn định, sử dụng k > n*(n+3)/2

# Lấy một điểm gần tâm khối lập phương
điểm = 0,5 + rand(n)*1e-3;
dữ liệu = rand(n, 10^4);
kdtree = KDTree(dữ liệu);
idxs, dists = knn(kdtree, point, k, true);

# Coords của k-Hàng xóm gần nhất
X = dữ liệu[:,idxs];

# Công thức bình phương nhỏ nhất cho hệ số
 C = điểm * cái(1,k); # nút trung tâm
dX = X - C # khác với nút trung tâm
 G = dX' * dX;
 F = G .* G;
 v = đường chéo(G);
 N = pinv(G) * G;
 N = mắt(N) - N;
 a = N * pinv(F*N) * v; đây là các hệ số

# Sử dụng phân bố nhiệt độ T = 25,4 * r^2
# có gradient phân tích là gradT = 25,4 * 2*x
X2 = X .* X;
C2 = C .* C;
T = 25,4 * n * trung bình(X2, 1)';
Tc = 25,4 * n * Mean(C2, 1)';
dT = T - Tc # khác với nút trung tâm

y = dX * (a .* dT); # Độ dốc ước tính
g = 2 * 25,4 * điểm # Phân tích

# in kết quả
@printf "Mức độ ước tính = %s\n" chuỗi(y')
@printf "Cấp độ phân tích = %s\n" chuỗi(g')
@printf "Lỗi tương đối = %.8f\n" vecnorm(gy)/vecnorm(g)

Phương pháp này có sai số tương đối khoảng 1%. Đây là kết quả của một số lần chạy...

Mức độ ước tính = [25,51670916224472 25,421038632006926 25,6711949674633]
Cấp độ phân tích = [25,41499027802736 25,44913042322385 25,448202594123806]
Lỗi tương đối = 0,00559934

Mức độ ước tính = [25,310574056859014 25,549736360607493 25,368056350800604]
Cấp độ phân tích = [25,43200914200516 25,43243178887198 25,45061497749628]
Lỗi tương đối = 0,00426558

gia hạn
Tôi không biết rõ về Python nhưng đây là bản dịch có vẻ hiệu quả

nhập numpy dưới dạng np
từ nhập scipy.spatial KDTree

n = 3;
k = 32;

# lấp đầy khối lập phương bằng các điểm ngẫu nhiên
dữ liệu = np.random.Rand(10000,n)
kdtree = KDTree(dữ liệu)

# chọn một điểm (ở tâm khối lập phương)
điểm = 0,5 * np.ones((1,n))

#Coords của k-Hàng xóm gần nhất
dists, idxs = kdtree.query(point, k)
idxs = idxs[0]
X = dữ liệu[idxs,:]

#Tính hệ số
C = (np.dot(point.T, np.ones((1,k)))).T # nút trung tâm
dX= X - C # khác với nút trung tâm
G = np.dot(dX, dX.T)
F = np.multiply(G, G)
v = np.diag(G);
N = np.dot(np.linalg.pinv(G), G)
N = np.eye(k) - N;
a = np.dot(np.dot(N, np.linalg.pinv(np.dot(F,N))), v) # đây là các hệ số

# Phân bố nhiệt độ là T = 25,4 * r^2
X2 = np.multiply(X, X)
C2 = np.multiply(C, C)
T = 25,4 * n * np.mean(X2, 1).T
Tc = 25,4 * n * np.mean(C2, 1).T # nút trung tâm
dT = T - Tc # khác với nút trung tâm

# Độ dốc phân tích ==> gradT = 2*25.4* x
g = 2 * 25,4 * điểm;
print( "g[]: %s" % (g) )

#Độ dốc ước tính
y = np.dot(dX.T, np.multiply(a, dT))
print( "y[]: %s, Lỗi tương đối = %.8f" % (y, np.linalg.norm(gy)/np.linalg.norm(g)) )

Cập nhật số 2
Tôi nghĩ tôi có thể sử dụng ASCII được định dạng thay vì LaTeX để viết thứ gì đó dễ hiểu hơn...

`Cho một tập hợp M vectơ có n chiều (gọi chúng là b_k), hãy tìm một tập hợp` hệ số (gọi chúng là a_k) mang lại ước tính tốt nhất về ma trận danh tính và vectơ 0`` M` (1) phút ||E - I||, trong đó E = tổng a_k b_k b_k` a_k k=1`` M` (2) phút ||z - 0||, trong đó z = tổng a_k b_k` a_k k=1```Lưu ý rằng các vectơ cơ sở {b_k} không bắt buộc`phải chuẩn hóa, trực giao hoặc thậm chí độc lập tuyến tính.``Đầu tiên, xác định các đại lượng sau:`` B ==> ma trận có các cột là b_k` G = B'.B ==> hoán vị của B nhân B` F = G @ G ==> @ đại diện cho sản phẩm hadamard` v = diag(G) ==> vectơ bao gồm các phần tử diag của G` `Điều trên các phép giảm thiểu tương đương với bài toán ràng buộc tuyến tính này Giải Fa = v` st Ga = 0 Đặt {X} biểu thị nghịch đảo Moore-Penrose của X.`Sau đó, nghiệm của bài toán tuyến tính có thể được viết:` N = I - {G} .G ==> chiếu vào không gian rỗng của G` a = N . {FN} . Tiện ích của các hệ số này là chúng cho phép bạn viết`các biểu thức rất đơn giản cho đạo hàm của một tensor field.```Cho D là toán tử del (hoặc nabla)` và d là toán tử sai phân viết nút trung tâm (hay còn gọi là số 0),` do đó, với bất kỳ đại lượng vô hướng/vectơ/tensor Y nào, chúng ta có:` dY = Y - Y_0``Cho x_k là vectơ vị trí của nút thứ k.`Và đối với các vectơ cơ sở của chúng ta, take` b_k = dx_k = x_k - x_0.``Giả sử rằng mỗi nút có một giá trị trường được liên kết với nó` (ví dụ: nhiệt độ) và giả sử mô hình bậc hai [khoảng x = x_0]` cho trường [g=gradient, H=hessian, :// là tích hai dấu chấm]`` Y = Y_0 + (x-x_0).g + (x-x_0)(x-x_0):H/2` dY = dx.g + dxdx:H/2` D2Y = I:H ==> Laplacian của Y```Đánh giá mô hình tại nút thứ k `` dY_k = dx_k.g + dx_k dx_k:H/2``Nhân với a_k và sum`` MMM` sum a_k dY_k = sum a_k dx_k.g + tổng a_k dx_k dx_k:H/2` k=1 k=1 k=1`` = 0.g + I:H/2` = D2Y / 2``Do đó, chúng ta có ước tính bậc hai của Laplacian`` M` Lap(Y_0) = tổng 2 a_k dY_k` k=1```Bây giờ hãy chơi cùng một trò chơi với mô hình tuyến tính` dY_k = dx_k.g``Nhưng lần này nhân với (a_k dx_k) và sum`` MM` sum a_k dx_k dY_k = sum a_k dx_k dx_k.g` k=1 k=1`` = Ig` = g```Nói chung, đạo hàm tại nút trung tâm có thể được ước tính là`` M` D#Y = sum a_k dx_k# dY_k` k=1`` M` D2Y = sum 2 a_k dY_k` k=1`` trong đó` # là viết tắt của tích {dot, cross, or tensor}` mang lại {div,curl, or grad} của Y` và` D2Y là viết tắt của Laplacian của Y` D2Y = D.DY = Lap(Y)

Về python - tính toán độ dốc 3D với khoảng cách điểm không đồng đều, chúng tôi đã tìm thấy một câu hỏi tương tự trên Stack Overflow: https://stackoverflow.com/questions/40338386/