有没有一种为 scipy.optimize.root 或 scipy.optimize.fsolve 动态创建非线性方程的好方法？答案

【问题标题】：Is there a good way to dynamically create non-linear equations for scipy.optimize.root or scipy.optimize.fsolve?有没有一种为 scipy.optimize.root 或 scipy.optimize.fsolve 动态创建非线性方程的好方法？
【发布时间】：2019-03-14 23:54:10
【问题描述】：

我需要求解一个大型非线性方程组（静态桁架系统）。这些方程源自节点 (xyz) 及其约束（位置、力）。

目前我们正在使用 Mathematica 来完成这项任务，但我们想迁移到 Python。但是使用 Mathematica（或 EES（工程方程求解器）或 SymPy）非常方便。我将一堆东西放在节点位置或节点上的力中，它会发挥一些作用，并自行创建方程式并结合输入并求解它们。

如果我想使用scipy.optimize.root，我必须以某种方式获得方程式。

scipy.optimize.root 和scipy.optimize.fsolve 需要以下格式的方程：

def func(x):
    out = [x[0]*cos(x[1]) - 4], 
           x[1]*x[0] - x[1] - 5)
    return out

但在我的例子中，将有多达 5000 个方程来定义系统。

我想到的一件事是使用eval() 并以某种方式将方程式摆弄成一个字符串。

最后，我想要一种面向对象的方法，其中节点或约束知道如何将自己转换为方程。一个非常简单的骨架可以是

n = Node(5, 2, 6)
n.to_equation()

f = ForceConstraint(1, 2, 3)
f.to_equation()

这会以某种方式转换为等式

x[0] - 5, 
x[1] - 2,
x[2] - 6,
x[2] ** 2 - x[1] * x[0] # and some non-linear stuff

描述整个系统。

基本上应该有一些神奇的部分来查看方程式和约束的匹配部分。例如。查看您在 Node1 的 x 方向上的所有信息并将其合并到方程中，或者搜索您在 Node2 上的 y 方向上的所有力的所有信息。

scipy 是完成这项工作的正确工具吗？有人知道如何做到这一点吗？

【问题讨论】：

方程式有多相似？它们有共同的结构吗？
不，它们不相似，如果节点相互依赖，可以只依赖一两个变量或数百个变量。
我看不出x[2] ** 2 - x[1] * x[0] 是如何从Node(5, 2, 6) 派生而来的。规则是什么？我认为您需要更详细地解释to_equation() 部分的外观，即更好地解释从输入到输出的转换。
我在问题中添加了一些信息。这个等式更像是一个例子。
您可以将 f 定义为 f(x, a, b, c,...)，并使用 scipy.optimize.root 中的 args=(a,b,c...) 关键字将 a、b、c 作为参数传递。参数可以是任何东西（布尔值、数组、数字）。也许你可以用它来思考一组不同的方程或参数？

标签： python python-3.x numpy scipy sympy

【解决方案1】：

我想你可能对symfit 感兴趣。这是我写的用来连接scipy 和sympy 的包。

我不确定您的具体方程式是什么，但您可以在sympy 中编写的任何表达式原则上都可以提供给symfit。例如，对于上面的简单示例，您可以编写：

from symfit import parameters, variables, Fit
import numpy as np

x0, x1, x2 = parameters('x0, x1, x2')
y0, y1, y2, y3 = variables('y0, y1, y2, y3')

model_dict = {
    y0: x0 - 5,
    y1: x1 - 2,
    y2: x2 - 6,
    y3: x2 ** 2 - x1 * x0
}

fit = Fit(model_dict, y0=np.array(0.0), y1=np.array(0.0), y2=np.array(0.0), y3=np.array(0.0))
fit_result = fit.execute()
print(fit_result)

symfit 中的Variable 和Parameter 对象只是sympy Symbol 的子类，因此您可以对这些表达式进行您想要的所有sympy 操作。例如，您可以将您的节点定义为

>>> x, x_0 = symbols('x, x_0')
>>> Node = x - x_0

然后通过重复应用来制作模型的线条，例如

>>> Node.subs({x: x1, x_0: 2})
x1 - 2

最后你添加你的约束和presto：适合模型！查看docs 了解更多信息或询问我任何后续问题。

【讨论】：