背景。
我正在嘗試寫一個 Python在Math SE上實作這個答案。您可能會發現以下背景很有用。
問題
我有一個實驗裝置，由三 (3) 個具有已知位置的接收器[xi, yi, zi]和一個具有未知位置的[x,y,z]發射器組成，該發射器以已知速度發射信號v。該信號在已知時間到達接收器ti。發射時間 , t,未知。
我希望找到到達角（即發射機的極坐標theta和phi），僅給出此資訊。
解決方案
它是不是能夠找到確切發射只有三個（3）接收器，除了獨特的情況下，少數（有跨越幾個偉大的答案數學SE解釋為什么是這種情況）。通常，至少需要四個（實際上是 >>4 個）接收器來唯一地確定發射器的直角坐標。
然而，可以“可靠地”估計到發射機的方向。設vi是表示接收器位置的向量i，ti是信號到達接收器的時間i，n是表示指向發射器（近似）方向的單位向量的向量，我們得到以下等式：
<n, vj - vi> = v(ti - tj)
（其中< >表示標量積）
...對于所有索引對i, j. 與 一起|n| = 1，系統通常有 2 個解決方案，通過平面反射對稱vi/vj/vk。然后我們可以通過簡單地在極坐標中寫入來確定phi和。thetan
執行。
我試圖寫一個 Python 實作上述解決方案，使用 scipy的fsolve。
from dataclasses import dataclass
import scipy.optimize
import random
import math
c = 299792
@dataclass
class Vertexer:
roc: list
def fun(self, var, dat):
(x,y,z) = var
eqn_0 = (x * (self.roc[0][0] - self.roc[1][0])) (y * (self.roc[0][1] - self.roc[1][1])) (z * (self.roc[0][2] - self.roc[1][2])) - c * (dat[1] - dat[0])
eqn_1 = (x * (self.roc[0][0] - self.roc[2][0])) (y * (self.roc[0][1] - self.roc[2][1])) (z * (self.roc[0][2] - self.roc[2][2])) - c * (dat[2] - dat[0])
eqn_2 = (x * (self.roc[1][0] - self.roc[2][0])) (y * (self.roc[1][1] - self.roc[2][1])) (z * (self.roc[1][2] - self.roc[2][2])) - c * (dat[2] - dat[1])
norm = math.sqrt(x**2 y**2 z**2) - 1
return [eqn_0, eqn_1, eqn_2, norm]
def find(self, dat):
result = scipy.optimize.fsolve(self.fun, (0,0,0), args=dat)
print('Solution ', result)
# Crude code to simulate a source, receivers at random locations
x0 = random.randrange(0,50); y0 = random.randrange(0,50); z0 = random.randrange(0,50)
x1 = random.randrange(0,50); x2 = random.randrange(0,50); x3 = random.randrange(0,50);
y1 = random.randrange(0,50); y2 = random.randrange(0,50); y3 = random.randrange(0,50);
z1 = random.randrange(0,50); z2 = random.randrange(0,50); z3 = random.randrange(0,50);
t1 = math.sqrt((x0-x1)**2 (y0-y1)**2 (z0-z1)**2)/c
t2 = math.sqrt((x0-x2)**2 (y0-y2)**2 (z0-z2)**2)/c
t3 = math.sqrt((x0-x3)**2 (y0-y3)**2 (z0-z3)**2)/c
print('Actual coordinates ', x0,y0,z0)
myVertexer = Vertexer([[x1,y1,z1], [x2,y2,z2], [x3,y3,z3]])
myVertexer.find([t1,t2,t3])
不幸的是，我在C/C 使用GSL.scipy之類的。我收到錯誤：
TypeError: fsolve: there is a mismatch between the input and output shape of the 'func' argument 'fun'.Shape should be (3,) but it is (4,).
...這似乎表明fsolve需要一個方形系統。
我該如何解決這個矩形系統？我似乎找不到任何有用的東西scipy 檔案。
如有必要，我愿意使用其他（Python）庫。
uj5u.com熱心網友回復：
正如您已經提到的，fsolve期望一個具有 N 個變數和 N 個方程的系統，即它找到函式 F 的根：R^N -> R^N。由于您有四個方程，您只需添加第四個變數。另請注意，這fsolve是一個遺留功能，建議改用它root。最后但并非最不重要的一點sqrt(x^2 y^2 z^2) = 1是，請注意，x^2 y^2 z^2=1當逼近 F 的雅可比時，后者更不容易受到由有限差分引起的舍入誤差的影響。
長話短說，你的班級應該是這樣的：
from scipy.optimize import root
@dataclass
class Vertexer:
roc: list
def fun(self, var, dat):
x,y,z, *_ = var
eqn_0 = (x * (self.roc[0][0] - self.roc[1][0])) (y * (self.roc[0][1] - self.roc[1][1])) (z * (self.roc[0][2] - self.roc[1][2])) - c * (dat[1] - dat[0])
eqn_1 = (x * (self.roc[0][0] - self.roc[2][0])) (y * (self.roc[0][1] - self.roc[2][1])) (z * (self.roc[0][2] - self.roc[2][2])) - c * (dat[2] - dat[0])
eqn_2 = (x * (self.roc[1][0] - self.roc[2][0])) (y * (self.roc[1][1] - self.roc[2][1])) (z * (self.roc[1][2] - self.roc[2][2])) - c * (dat[2] - dat[1])
norm = x**2 y**2 z**2 - 1
return [eqn_0, eqn_1, eqn_2, norm]
def find(self, dat):
result = root(self.fun, (0,0,0,0), args=dat)
if result.success:
print('Solution ', result.x[:3])
轉載請註明出處，本文鏈接：https://www.uj5u.com/ruanti/406022.html
標籤：
上一篇：找出兩個方程的可能值