用SciPy求解这个矩形的非线性系统

本教程将介绍用SciPy求解这个矩形的非线性系统的处理方法，这篇教程是从别的地方看到的，然后加了一些国外程序员的疑问与解答，希望能对你有所帮助，好了，下面开始学习吧。

问题描述

背景。

我正在尝试在上编写python答案的python实现。您可能会发现以下背景知识很有用。

我有一个实验设置，由三(3)个接收器和一个发射器组成，其中三(3)个接收器具有已知位置[xi, yi, zi]，位置[x,y,z]以速度v发射信号。该信号在已知时间ti到达接收器。发射时间t。

我只想找出到达角(即发射机的极坐标theta和phi)，仅给出此信息。

仅用三(3)个接收器中有几个很好的答案解释了为什么会这样)。通常，至少需要四个(实际上是&gt；&gt；4)接收器才能唯一确定发射器的直角坐标。

但是，可以可靠地估计到发射机的。设vi是表示接收器位置的矢量i，ti是信号到达接收器的时间i，n是表示指向发射机(近似)方向的单位矢量的矢量，我们得到以下方程：

<n, vj - vi> = v(ti - tj)

对于所有索引对i、j。与|n| = 1一起，系统一般有两个解，通过vi/vj/vk在平面上反射对称。然后，我们只需在极坐标中书写n即可确定phi和theta。

实现。

我已尝试使用scipy的fsolve编写上述解决方案的python实现。

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++中解决此类问题的经验要丰富得多，而使用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)库。

推荐答案

正如您已经提到的，fsolve期望一个具有N个变量和N个方程的系统，即它找到函数F：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])

好了关于用SciPy求解这个矩形的非线性系统的教程就到这里就结束了，希望趣模板源码网找到的这篇技术文章能帮助到大家，更多技术教程可以在站内搜索。