计算求解 100 个变量的非线性方程组答案

【问题标题】：Solving System of Nonlinear Equations in 100 variables Computationally计算求解 100 个变量的非线性方程组
【发布时间】：2021-10-03 17:37:07
【问题描述】：

我一直在尝试使用 Python 中的 Sympy/Numpy 求解一个由 10,000 个非线性方程组组成的系统，该方程组包含 100 个变量。

目前的尝试：

我在 Sympy 中尝试了 nsolve 和 solve，在 numpy.linalg 中求解，所有这些都在等待 5-6 小时后仍在运行（我最终在 RAM 用完时强制停止它们）。

使用 Sympy 本身生成一组方程大约需要 1 小时。我切换到 SageMath（Windows 原生），它似乎可以更好地生成方程本身（约 3 分钟），但仍然无法解决它们。

有没有办法使用 SageMath/Python 本身中的任何特定语言或技巧来优化运行，或者我应该寻找更强大的系统来运行代码？

我的系统是 i7-11300H/16GB RAM。

编辑：linalg.solve 是一个错误，因为我最初认为它是一个线性系统，后来意识到它不是。

编辑：

from sympy import *
x,t=symbols('x,t')
a11, a12, a13, a14, a15, a16, a17, a18, a19, a21, a22, a23, a24 = symbols('a11, a12, a13, a14, a15, a16, a17, a18, a19, a21, a22, a23, a24')
Coeffs = [a11, a12, a13, a14, a15, a16, a17, a18, a19, a21, a22, a23, a24]
#E = expression in x,t,aij, 0<i<11,0<j<11 (nonlinear in aij)
## Sample
E = 1/2*(9*(8*t + 1)*a11 + 3*(t**2 + 2*t + 1)*a21 + 4*(t**3 + 3*t**2 + 3*t + 1)*a12 + 5*(t**4 + 4*t**3)*a13 + 6*(t**5 + 5*t**4)*a14 + 7*(t**6 + 6*t**5 + 1)*a15 + 8*(t**7 + 7*t)*a16 + 2*a17*(t + 1) + a18)*x**2 + 1/48*(13860*(252*(t**10 + 10*t)*a19 + 1260*(t**2 + 2*t)*a21 + 840*(t**3 + 3*t**2 + 3*t)*a22 + 630*(t**4)*a23 + 504*(t**5 + 10*t**2 + 5*t))*a24)
Eqs = [E.subs({x:1/k,t:1/m}) for k in range(1,100) for m in range(1,100)]
sol = solve(Eqs, Coeffs)

还尝试使用 0 数组作为初始值的 nsolve。

【问题讨论】：

对于为什么这些事情运行缓慢的原因有很多可能的解释，并且可能有很多方法可以让您进行更改以加快速度。你没有给我足够的细节让我说更具体的事情。我也不清楚，因为您说系统是非线性的，但您使用的是 linalg.solve ，它用于线性方程组的平方系统。最后，令我惊讶的是，您期望从这种过度约束的系统中找到任何解决方案。我从来没有见过 linalg.solve 只要你描述的那样接近任何地方，但我想我不明白你在做什么。
@OscarBenjamin linalg.solve 是一个错误，是的。我最初将它们视为线性系统，后来意识到我正在处理非线性，切换回 sympy。
另外，要么你没有任何解，要么你有许多可以简化的共线方程。而且我看不出 a_ij 变量的非线性在哪里。
确实，系统看起来非常冗余。我的猜测是您正试图以迂回的方式将x,t 项的系数等同于零。我添加了一个答案，显示如何直接执行此操作。

标签： python numpy math sympy sage

【解决方案1】：

在你更新你的问题后，我想我明白你在做什么，但我认为你没有以正确的方式解决问题。

首先，您定义方程组的方式引入了浮点系数，并且带有浮点数的多项式方程可能非常病态，因此请确保使用例如S(1)/2 或 Rational(1, 2) 而不是 1/2 像这样：

from sympy import *
x,t=symbols('x,t')
Coeffs = symbols('a11, a12, a13, a14, a15, a16, a17, a18, a19, a21, a22, a23, a24')
a11, a12, a13, a14, a15, a16, a17, a18, a19, a21, a22, a23, a24 = Coeffs

#E = expression in x,t,aij, 0<i<11,0<j<11 (nonlinear in aij)
## Sample
E = (
      S(1)/2*(9*(8*t + 1)*a11
    + 3*(t**2 + 2*t + 1)*a21
    + 4*(t**3 + 3*t**2 + 3*t + 1)*a12
    + 5*(t** 4 + 4*t**3)*a13
    + 6*(t**5 + 5*t**4)*a14
    + 7*(t**6 + 6*t**5 + 1)*a15
    + 8*(t**7 + 7*t)*a16
    + 2*a17*(t + 1) + a18)*x**2
    + S(1)/48*(13860*(
          252*(t**10 + 10*t)*a19
        + 1260*(t**2 + 2*t)*a21
        + 840*(t**3 + 3*t**2 + 3*t)*a22
        + 630*(t**4)*a23
        + 504*(t**5 + 10*t**2 + 5*t) )*a24
        )
    )

所以你有E，看起来像这样：

In [2]: E
Out[2]: 
    ⎛          ⎛     10         ⎞             ⎛      2         ⎞             ⎛     3         2     
a₂₄⋅⎝13860⋅a₁₉⋅⎝252⋅t   + 2520⋅t⎠ + 13860⋅a₂₁⋅⎝1260⋅t  + 2520⋅t⎠ + 13860⋅a₂₂⋅⎝840⋅t  + 2520⋅t  + 25
───────────────────────────────────────────────────────────────────────────────────────────────────
                                                                                48                 

    ⎞                4            5             2             ⎞      ⎛                     ⎛   3   
20⋅t⎠ + 8731800⋅a₂₃⋅t  + 6985440⋅t  + 69854400⋅t  + 34927200⋅t⎠    2 ⎜a₁₁⋅(72⋅t + 9)   a₁₂⋅⎝4⋅t  + 
─────────────────────────────────────────────────────────────── + x ⋅⎜────────────── + ────────────
                                                                     ⎝      2                      

    2           ⎞       ⎛   4       3⎞       ⎛   5       4⎞       ⎛   6       5    ⎞       ⎛   7   
12⋅t  + 12⋅t + 4⎠   a₁₃⋅⎝5⋅t  + 20⋅t ⎠   a₁₄⋅⎝6⋅t  + 30⋅t ⎠   a₁₅⋅⎝7⋅t  + 42⋅t  + 7⎠   a₁₆⋅⎝8⋅t  + 
───────────────── + ────────────────── + ────────────────── + ────────────────────── + ────────────
  2                         2                    2                      2                      2   

    ⎞                           ⎛   2          ⎞⎞
56⋅t⎠                 a₁₈   a₂₁⋅⎝3⋅t  + 6⋅t + 3⎠⎟
───── + a₁₇⋅(t + 1) + ─── + ────────────────────⎟
                       2             2          ⎠

E 和 a24 的第一项使这个非线性，但只是轻微的，因为它是二次和多项式，而不是超越或其他东西（你以前只将你的方程描述为非线性，这可能意味着任何东西）。

您正在为 x 和 t 中的每一个替换 100 个不同的值，然后尝试求解 10000 个方程以使 E 中的每个值都为零。几乎可以保证，这些方程的唯一可能解是那些使x,t 多项式的所有系数为零的解。我猜你实际上想要做的是找到使所有E 为零的ai 符号的值x,t 但我们不需要10000 个方程来做到这一点。我们可以提取系数并将其求解为方程：

In [3]: eqs = Poly(E, [x, t]).coeffs()

In [4]: eqs
Out[4]: 
⎡       7⋅a₁₅                  5⋅a₁₃                                   3⋅a₂₁                       
⎢4⋅a₁₆, ─────, 3⋅a₁₄ + 21⋅a₁₅, ───── + 15⋅a₁₄, 2⋅a₁₂ + 10⋅a₁₃, 6⋅a₁₂ + ─────, 36⋅a₁₁ + 6⋅a₁₂ + 28⋅a
⎣         2                      2                                       2                         

                  9⋅a₁₁           7⋅a₁₅         a₁₈   3⋅a₂₁                             363825⋅a₂₃⋅
₁₆ + a₁₇ + 3⋅a₂₁, ───── + 2⋅a₁₂ + ───── + a₁₇ + ─── + ─────, 72765⋅a₁₉⋅a₂₄, 145530⋅a₂₄, ───────────
                    2               2            2      2                                     2    

a₂₄                                                                                                
───, 242550⋅a₂₂⋅a₂₄, 363825⋅a₂₁⋅a₂₄ + 727650⋅a₂₂⋅a₂₄ + 1455300⋅a₂₄, 727650⋅a₁₉⋅a₂₄ + 727650⋅a₂₁⋅a₂₄
                                                                                                   

                              ⎤
 + 727650⋅a₂₂⋅a₂₄ + 727650⋅a₂₄⎥
                              ⎦

In [5]: %time [sol] = solve(eqs, Coeffs, dict=True)
CPU times: user 236 ms, sys: 8 ms, total: 244 ms
Wall time: 244 ms

In [6]: sol
Out[6]: 
⎧     a₁₈                                               -4⋅a₁₈                 ⎫
⎨a₁₁: ───, a₁₂: 0, a₁₃: 0, a₁₄: 0, a₁₅: 0, a₁₆: 0, a₁₇: ───────, a₂₁: 0, a₂₄: 0⎬
⎩      63                                                  7                   ⎭

In [7]: checksol(eqs, sol)
Out[7]: True

这表明方程组是欠定的，例如a18 可以任意提供 a11 = a18/63 等。如果解决方案字典中未列出其中一个未知数，则它可以独立于其他取任意值。

我在那里添加了一个%time，以表明在我不是很快的计算机上解决这个问题需要 244 毫秒（在当前的 sympy master 分支上——使用 sympy 1.8 大约需要 0.5 秒）。请注意，我确实安装了 gmpy2，它可以使 SymPy 中的各种事情变得更快。安装 sympy 1.9rc1 和 gmpy2（pip install --pre --upgrade sympy 和 pip install gmpy2）可能会加快你的速度。还值得注意的是，最近在 sympy 中修复了此类系统的一些错误，因此对于其他示例，1.9 可能会给出更准确的结果：

https://github.com/sympy/sympy/pull/21883

【讨论】：

如前所述，这只是一个简单的方程式。 Orignal 涉及所有 100 个 aijs。 eqs = Poly(E, [x,t]).coeffs() 与原始方程给出 [E]
大概那么您的实际E 不包括符号x 和t，您需要更改该行。如果您无法解释，我不会尝试猜测您在做什么。

【解决方案2】：

您可以将线程用于您的目的。

这里是 Corey Schäfer 的解释 --> https://youtu.be/IEEhzQoKtQU

核心思想是并行运行代码，这意味着您的代码不会像“......计算第一个代码块......（完成）......计算第二个代码块......完成”等等在。通过线程，您的代码将一次运行多个代码块。

一个好的 Python 库是 multiprocessing 库。

但首先要确定您的机器可以使用多少个处理器。

import multiprocessing as mp
print("Number of processors: ", mp.cpu_count())

--> Number of processors:  4

例如，我的电脑只有 4 个可用处理器

这里是 python 中的线程示例：

import requests
from multiprocessing.dummy import Pool as ThreadPool 

urls = ['https://www.python.org',
        'https://www.python.org/about/']

def get_status(url):
    r = requests.get(url)
    return r.status_code

if __name__ == "__main__":
    pool = ThreadPool(4)  # Make the Pool of workers
    results = pool.map(get_status, urls) #Open the urls in their own threads
    pool.close() #close the pool and wait for the work to finish 
    pool.join()

上面的代码将循环遍历 URL 列表并输出 URL 的状态代码。如果没有线程，代码将以迭代的方式执行此操作。第一个 URL --> 状态 200 ... NEXT ... 第二个 URL --> 状态 200 ...NEXT ... 等等。

希望能帮到你。

【讨论】：