在 Python 上求解非线性系统

Question

我通常依靠 Wolfram Mathematica 来做这类事情，但最近我一直在研究 Python，结果好多了。基本上，我正在为如下系统寻找数值解决方案。

系统：

嗯，我知道有解决方案，因为 Wolfram Mathematica 找到了一个 (0.0858875,0.0116077,-0.156661,1.15917)。我在 Python 中尝试做的是这个暴力代码。

import numpy as np

START = -3
END = 3
STEP = 0.1

for r0 in np.arange(START, END, STEP):
    for r1 in np.arange(START, END, STEP):
        for r2 in np.arange(START, END, STEP):
            for r3 in np.arange(START, END, STEP):
                eq0 = r0*r2+r1*r3
                eq1 = r0*r1+r1*r2+r2*r3+r0*r3
                eq2 = r0**2+r1**2+r2**2+r3**2-4*(r0+r1+r2+r3)**2

                if (eq0 == 0 and eq1 < 0 and eq2 < 0): 
                    print(r0, r1, r2, r3)

编辑： 我同意 -0.00001< eq0 < 0.00001 而不是 eq1 == 0

好吧，虽然在这种情况下没有找到解决方案，但对于我正在处理的其他系统，蛮力方法运行良好，特别是当方程和变量较少时。从四个变量开始，真的很难。

对不起，如果我要求太多。我是 Python 的新手，所以我也不知道这是否真的微不足道。也许 fsolve 会有用？我不确定它是否适用于不平等现象。此外，即使我遇到的系统只有等式，它们的变量总是比方程多，就像这个：系统 2：
.

所以'fsolve'不合适吧？

Answer 1

不是真正的答案，但您可以使用 product

来简化这个问题

from itertools import product
import numpy as np


START = -3
END = 3
STEP = 0.1

for r0, r1, r2, r3 in product(np.arange(START, END, STEP), repeat=4):
    print(r0, r1, r2, r3)

Answer 2

不确定您的问题是寻根问题还是约束最小化问题，但请查看 scipy.optimize 和 scipy.linprog，也许其中一种方法可以适应您的应用程序.

Answer 3

只要您的系统包含不等式，您就需要将其表述为优化问题，并使用 scipy.optimize.minimize. Otherwise, you can use scipy.optimize.root 或 scipy.optimize.fsolve 求解方程组。请注意，前者也正是 root 和 fsolve 中幕后所做的，即两者都解决了最小二乘优化问题。

一般来说，问题

g_1(x) = 0, ..., g_m(x) = 0
h_1(x) < 0, ..., h_p(x) < 0

可以表述为

min g_1(x)**2 + ... + g_m(x)**2

s.t. -1.0*(h_1(x) + eps) >= 0
                  .
                  .
     -1.0*(h_p(x) + eps) >= 0

其中 eps 是对严格不等式建模的容差。

因此，您可以按如下方式解决您的第一个问题：

import numpy as np
from scipy.optimize import minimize

def obj(r):
    return (r[0]*r[2]+r[1]*r[3])**2

eps = 1.0e-6

constrs = [
    {'type': 'ineq', 'fun': lambda r: -1.0*(r[0]*r[1] + r[1]*r[2] + r[2]*r[3] + r[0]*r[3] + eps)},
    {'type': 'ineq', 'fun': lambda r: -1.0*(np.sum(r**2) - 4*(np.sum(r))**2 + eps)}
]

# res.x contains the solution
res = minimize(obj, x0=np.ones(4), constraints=constrs)

你的第二个问题可以用类似的方法解决。在这里，您只需要删除约束即可。或者，您可以在值得一提的地方使用 root，它为函数 F: R^N -> R^N 求解 F(x) = 0，即 N 变量的函数 returns 和 N维向量。如果您的函数的方程式少于变量，您可以简单地用零填充向量：

import numpy as np
from scipy.optimize import root

def F(r):
    vals = np.zeros(r.size)
    vals[0] = np.dot(r[:5], r[1:]) + r[0]*r[5]
    vals[1] = r[0]*r[3] + r[1]*r[4] + r[2]*r[5]
    vals[2] = np.sum(r**2) - 3*np.sum(r)**2
    return vals

# res.x contains your solution
res = root(F, x0=np.ones(6))

在 Python 上求解非线性系统

Solving a non-linear system on Python

python

optimization

linear-algebra

nonlinear-optimization