python - 在Python中求解非线性方程: the answer is the same as initial guess

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) 看到这种情况的发生——它返回 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 用于解析求解方程。它不关心一堆 float 。给它一个符号方程:

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 收敛失败的预期:方程没有真正的解。