在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 在这个方程中遇到的问题是重要性的丧失。您将 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用于解析方程。它不关心一堆浮点数。给它一个符号方程：

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收敛失败的期望：方程没有真正的解决方案。