scipy Python中的三变量指数方程-求解器选项

mspsb9vt 于 2023-08-05 发布在 Python

关注(0)|答案(2)|浏览(149)

问这里是我最后的手段，我不知道还能在哪里寻找或我的下一个选择将是什么，除了改变语言。
我有一组3重指数方程，有3个变量。指数内以及指数外都有未知变量。
我试过sympy的solve、nsolve和nonlinsolve以及scipy的fsolve。Scipy的fsolve给出了一个解，但是这个解包含负值，这是不允许的，因为我试图计算的是一个人生病的过渡强度，它不可能是负值。我的猜测是基于经验数据，应用任何其他随机猜测没有一个健全的理由背后的决定。Sympy的solve返回一个[]或只是继续运行，这取决于是从main还是独立调用（稍后复制），nsolve会吐出下面的错误，而nonlinsolve会无休止地运行而不返回任何东西。
x1c 0d1x的数据
我在主题列表中选择了MATLAB作为工具，但是，我宁愿留在Python中，因为我对MATLAB非常生疏，并且我已经写了一堆代码（无论多么糟糕，因为我不是一个真实的的程序员，只是一个使用Python写论文的人）来达到这个不愉快的点，我现在被卡住了。
如果有帮助的话，这就是方程的内容，以及我正在进行的求解调用的示例（从代码的其他部分提取）：

import sympy as sm
from sympy import nsolve
from sympy import solve

sigma_R,sigma_SU,sigma_MU=sm.symbols("sigma_R,sigma_SU,sigma_MU", real=True)

f1 = -0.262617004746481*sigma_SU*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/((-sigma_MU -   sigma_R - sigma_SU + 0.000931334226877499)*(0.260884896827677*sigma_MU*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00170434902663555) + 0.256617389979587*sigma_R*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00363094118807797) + 0.262617004746481*sigma_SU*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.000931334226877499) + 0.719582172729324*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU))) + 1.54

f3 = -0.260884896827677*sigma_MU*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/((-sigma_MU - sigma_R - sigma_SU + 0.00170434902663555)*(0.260884896827677*sigma_MU*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00170434902663555) + 0.256617389979587*sigma_R*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00363094118807797) + 0.262617004746481*sigma_SU*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.000931334226877499) + 0.719582172729324*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU))) + 0.3

f2 = -0.256617389979587*sigma_R*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/((-sigma_MU - sigma_R - sigma_SU + 0.00363094118807797)*(0.260884896827677*sigma_MU*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00170434902663555) + 0.256617389979587*sigma_R*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00363094118807797) + 0.262617004746481*sigma_SU*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.000931334226877499) + 0.719582172729324*sm.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU))) + 0.37

sol=solve((f1,f2,f3),(sigma_R,sigma_SU,sigma_MU),simplify=False)
print(sol)

字符串
我的要求是：
1.有没有一种方法我可以调试，步骤等。在任何solve函数中，它是一个黑盒的事实并没有帮助，因为如果我以某种方式发送错误的数据，我不会得到任何例外。
1.如果有人知道我没有尝试过的方法，请分享;我没办法了。即使是一个数学建议也会起作用，但我看不出有任何方法可以简化方程本身。
如上所述，我尝试过：

sol = solve((f1,f2,f3),(sigma_R,sigma_SU,sigma_MU),simplify=False)
sol = nsolve((f1, f2,f3), (sigma_R,sigma_SU,sigma_MU), (0.0000813,0.0000202,0.000384))
sol = sm.nonlinsolve([y0,y1,y2],[sigma_R,sigma_SU,sigma_MU])

型
fsolve的部分：

from scipy.optimize import fsolve
import numpy as np
import math
import datetime as dt

def equations(z):
    sigma_SU=z[0]
    sigma_MU=z[1]
    sigma_R=z[2]
    f=np.empty(3)
    f[0] = -0.262617004746481*sigma_SU*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/((-sigma_MU - sigma_R - sigma_SU + 0.000931334226877499)*(0.260884896827677*sigma_MU*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00170434902663555) + 0.256617389979587*sigma_R*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00363094118807797) + 0.262617004746481*sigma_SU*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.000931334226877499) + 0.719582172729324*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU))) + 1.54
    f[1] = -0.260884896827677*sigma_MU*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/((-sigma_MU - sigma_R - sigma_SU + 0.00170434902663555)*(0.260884896827677*sigma_MU*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00170434902663555) + 0.256617389979587*sigma_R*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00363094118807797) + 0.262617004746481*sigma_SU*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.000931334226877499) + 0.719582172729324*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU))) + 0.3
    f[2] =  -0.256617389979587*sigma_R*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/((-sigma_MU - sigma_R - sigma_SU + 0.00363094118807797)*(0.260884896827677*sigma_MU*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00170434902663555) + 0.256617389979587*sigma_R*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.00363094118807797) + 0.262617004746481*sigma_SU*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU)/(-sigma_MU - sigma_R - sigma_SU + 0.000931334226877499) + 0.719582172729324*np.exp(-45*sigma_MU - 45*sigma_R - 45*sigma_SU))) + 0.37
    return f

myGuess=np.array([0.0000813,0.0000202,0.000384])
z=np.array([0,0,0])
counter=0

ct1 = dt.datetime.now()
print(ct1)

z=fsolve(equations, myGuess)

ct2=dt.datetime.now()
print(ct2-ct1)
print(z)
print(counter)

型

scipy

来源：https://stackoverflow.com/questions/76699766/three-variables-exponential-equations-in-python-solver-options

2条答案

按热度按时间

i2loujxw1#

你需要重写你的公式，以确保合理性和性能，从而正确使用Numpy。然后，从fsolve（MINPACK hybr*）切换到minimize。初始向量（0，0，0）似乎比您提供的猜测更容易收敛。在以下四种无需额外“帮助”即可成功收敛的方法中，

Nelder-Mead
鲍威尔
L-BFGS-B，和
trust-constr，

trust-constr是最快收敛到与机器精度“精确”的解的方法：

from scipy.optimize import minimize
import numpy as np

c1 = np.array((0.262617004746481000, 0.260884896827677000, 0.256617389979587000))
c2 = np.array((0.000931334226877499, 0.001704349026635550, 0.003630941188077970))
c3 = np.array((1.54, 0.30, 0.37))

def equations(z: np.ndarray) -> np.ndarray:
    # sigma_SU, sigma_MU, sigma_R = z

    inner1 = -z.sum()
    inner2 = np.exp(45*inner1)
    inner3 = c1 * z * inner2/(inner1 + c2)
    f = c3 - inner3/(inner3.sum() + 0.719582172729324 * inner2)
    return f

def lstsq_cost(z: np.ndarray) -> float:
    f = equations(z)
    return f.dot(f)

def regression_tests() -> None:
    expected = (1.48037098, 0.29461646, 0.33099203)
    actual = equations(np.array((8.13e-5, 2.02e-5, 3.84e-4)))
    assert np.allclose(expected, actual)

    expected = (0.93253452, 0.21240985, 0.06805596)
    actual = equations(np.array((-0.49717498, -0.00991807,  0.50892057)))
    assert np.allclose(expected, actual)

def optimize() -> None:
    result = minimize(
        fun=lstsq_cost, method='powell',
        x0=(0, 0, 0), tol=1e-9,
    )
    assert result.success, result.message
    np.set_printoptions(precision=2)
    print('Function at roots:', equations(result.x))
    np.set_printoptions(precision=20)
    print('Roots:', result.x)

if __name__ == '__main__':
    regression_tests()
    optimize()

个字符
暴力外部搜索循环：

for method in (
        'Nelder-Mead',
        'Powell',
        'CG',
        'BFGS',
        'Newton-CG',
        'L-BFGS-B',
        'TNC',
        'COBYLA',
        'SLSQP',
        'trust-constr',
        'dogleg',
        'trust-ncg',
        'trust-exact',
        'trust-krylov',
    ):
        try:
            result = minimize(
                fun=lstsq_cost, method=method,
                x0=(0, 0, 0), tol=1e-9,
            )
            assert result.success, result.message
            y = equations(result.x)
            if y.dot(y) > 0.1:
                continue
            np.set_printoptions(precision=2)
            print(method, 'converged in nfev =', result.nfev)
            print('Function at roots:', y)
            np.set_printoptions(precision=20)
            print('Roots:', result.x)
            print()
        except Exception:
            pass

Nelder-Mead converged in nfev = 285
Function at roots: [-9.10e-08 -1.46e-07  3.42e-08]
Roots: [ 0.0032675556912995915   0.00011212905180447633 -0.001511370899412161  ]

Powell converged in nfev = 1135
Function at roots: [ 6.51e-14 -7.53e-14  5.76e-14]
Roots: [ 0.0032675561103559787   0.00011212902858510153 -0.0015113712768471598 ]

L-BFGS-B converged in nfev = 188
Function at roots: [2.39e-05 2.28e-05 7.52e-06]
Roots: [ 0.0032676040558970083   0.00011211616767734874 -0.0015114200558563654 ]

trust-constr converged in nfev = 292
Function at roots: [-0. -0.  0.]
Roots: [ 0.0032663289158950093   0.00011197597802975169 -0.0015103347529418544 ]

的字符串

赞(0）回复(0）举报 2023-08-05

yacmzcpb2#

由于我在评论中提到了它，为了完整性，我还将为@Reinderien的答案添加一个替代方法（尽管他们的答案更好，更完整，所以它应该比我的答案更好）。
如果你不想手动重写方程并得到等价的结果，你可以使用sympy的lambdify函数来为你正在使用的三个函数中的每一个生成向量化函数（这些函数在运行时不会是最有效的，但它是获得向量化函数的快速和肮脏的方法）。您可以使用“lambdified”函数来创建向量化的equations函数。可以使用scipy.optimize.root（fsolve是传统语法;你应该使用这个函数代替）与np.zeros(3)的初始猜测（在这种情况下，root不收敛于问题的初始猜测）。

import sympy as sp
import numpy as np
from scipy.optimize import root

sigma_R, sigma_SU, sigma_MU = sp.symbols("sigma_R,sigma_SU,sigma_MU", real=True)
variables = (sigma_SU, sigma_MU, sigma_R)

# define f1, f2, and f3 using sympy as in the question

f1_f = sp.lambdify((variables), f1)
f2_f = sp.lambdify((variables), f2)
f3_f = sp.lambdify((variables), f3)

def equations(z):
    return np.array([f1_f(*z), f2_f(*z), f3_f(*z)])

sol = root(equations, np.zeros(3))
print(sol)
print(sol.x)

字符串
输出量：

message: The solution converged.
success: True
 status: 1
    fun: [-6.737e-12 -2.788e-11  9.857e-12]
      x: [ 3.268e-03  1.121e-04 -1.511e-03]
   nfev: 51
   fjac: [[-8.232e-01  2.522e-01 -5.087e-01]
          [ 2.705e-01 -6.136e-01 -7.418e-01]
          [ 4.992e-01  7.483e-01 -4.369e-01]]
      r: [ 2.049e+03 -2.050e+03  2.379e+03  1.686e+03  3.289e+01
          -4.480e+02]
    qtf: [-2.165e-13  2.192e-10 -1.184e-09]

array([ 0.0032675561102804376 ,  0.00011212902858868879,
       -0.0015113712767817957 ])

型
P.S.如果你通过计算前后时间来计算函数的速度，请使用time库中的time.perf_counter()，而不是你使用的datetime方法。

from time import perf_counter

start = perf_counter()
# code here
end = perf_counter()
print(end - start)

型
如果你真的想认真对待时间，你可以看看timeit库。

赞(0）回复(0）举报 2023-08-05

我来回答

scipy Python中的三变量指数方程-求解器选项

2条答案

相关问题

热门标签

最新问答