求解python中的非线性方程：答案与初始猜测相同

Question

所以我有这个复杂的方程式需要求解。我认为最后 x 应该是 1E22 的顺序。但是这段代码的问题是它会导致我的整个系统崩溃。有解决办法吗？我尝试了 scipy.optimize.root 但它并没有真正解决这个数量级的任何问题（它给出了最终答案作为初始猜测而没有任何迭代）。

from scipy.optimize import fsolve
import math
import mpmath
import scipy
import sympy
from sympy.solvers import solve
from sympy import Symbol
from sympy import sqrt,exp

x = Symbol('x',positive=True)

cs = 507.643E-12
esi = 1.05E-10 
q = 1.6E-19
T = 300
k = 1.381E-23
ni = 1.45E16

print(solve(exp(x/((2*cs/(esi*q))**2)) - ((x/ni)**(esi*k*T)),x))

def func(N):
    return (math.exp(N/math.pow(2*cs/(esi*q),2)) - math.pow(N/ni,(esi*k*T)))

n_initial_guess = 1E21
n_solution = fsolve(func, n_initial_guess)

print ("The solution is n = %f" % n_solution)
print ("at which the value of the expression is %f" % func(n_solution))
print(scipy.optimize.root(func, 1E22,tol=1E-10))

scipy 函数均无效。 sympy 函数使我的笔记本电脑崩溃。 Matlab 是这个的理想选择吗？

Answer 1

SciPy

的数值解

SciPy 这个等式的问题是 loss of significance。您将 N 提高到 esi*k*T 的微小幂，这使得它非常接近 1；在浮点运算中，它恰好变为 1。类似地，来自指数的部分变为 1。然后两部分相减，留下 0 - 方程似乎已经解决了。您可以通过打印 func(1E21) 看到这种情况发生 - 它 returns 0.

处理失去意义的方法是重写方程，从原来的形式

exp(x/((2*cs/(esi*q))**2)) == (x/ni)**(esi*k*T)

通过提高双方的权力 1/(esi*k*T):

exp(x*esi*q**2/(2*cs*k*T)**2)) == x/ni

所以func变成了

def func(N):
    return np.exp(N*esi*q**2/(k*T*(2*cs)**2)) - (N/ni)

（建议将 NumPy 函数与 SciPy 求解器一起使用。）也就是说，求解器（例如 root(func, 1E10)）将报告无法收敛到解。

使用 SymPy 的符号解决方案

SymPy 用于分析求解方程。它不关心一堆浮点数。给它一个符号方程式：

x, a, b, c = symbols('x, a, b, c', positive=True)
sol = solve(exp(x/a) - (x/b)**c, x)[0]

得到的解为-c*LambertW(-a/(b*c))/a。那么就可以评价了：

cs = 507.643E-12
esi = 1.05E-10 
q = 1.6E-19
T = 300
k = 1.381E-23
ni = 1.45E16

print(sol.evalf(subs={a: (2*cs/(esi*q))**2, b: ni, c: esi*k*T}))

打印 -21301663061.0653 - 4649834682.69762*I 确认人们已经从与 SciPy 的收敛失败中预期的结果：方程没有真正的解。

求解python中的非线性方程：答案与初始猜测相同

Solving a non-linear equation in python: the answer is the same as initial guess

python

sympy

scipy

equation-solving

SciPy

使用 SymPy 的符号解决方案