百元非线性方程组的计算求解

Question

我一直在尝试使用Python中的Sympy/Numpy来解决一个包含10000个非线性方程组和100个变量的系统

迄今为止：

我尝试了Sympy的solve和solve，numpy.linalg中的solve和solve，经过5-6个小时的等待后，它们仍然在运行（当我的内存用完时，我强制停止了它们）

用Symphy本身生成方程组需要约1小时。我转向了SageMath（Windows native），它似乎能更好地生成方程（约3分钟），但仍然无法解决这些问题

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

我的系统是i7-11300H/16gbram

编辑：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

Answer 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          ⎠

带有a24的E的第一项使得这个非线性，但只是稍微非线性，因为它是二次和多项式，而不是说超越的或什么的（你之前只描述你的方程为非线性，可能意味着什么）

您将用100个不同的值替换x和t中的每一个，然后尝试求解10000个方程，使这些值的E为零。几乎可以保证，这些方程的唯一可能解是那些使x,t多项式的所有系数为零的解。我想你实际上想做的是找到ai符号的值，这些符号使E为零，但我们不需要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等。如果其中一个未知数没有列在解的dict中，那么它可以独立于其他未知数而取任意值

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

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

Answer 2

您可以根据需要使用线程

下面是Corey Schäfer的解释>https://youtu.be/IEEhzQoKtQU

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

一个好的python库是多处理库

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

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。。。下一个第二个URL>；状态200…下一个。。。等等

我希望我能帮你一点忙