在 python/SymPy 中数值求解不适定非线性方程答案

【问题标题】：Solving ill-posed non-linear equations numerically in python/SymPy在 python/SymPy 中数值求解不适定非线性方程
【发布时间】：2022-09-27 15:59:22
【问题描述】：

我试图通过运行下面的代码来获得解决方案。

Python 只是“挂起”，找不到数值解。我可以使用手机上的应用程序 (Desmos) 绘制函数图并轻松找到数值解，0.024。 python在解决小数点后2位时有限制吗？

import sympy

x = sympy.symbols(\'x\')
e_1 = x**-0.5
e_2 = -2*sympy.log(0.0001*3.7**-1*0.05**-1+2.51*350000**-1*x**-0.5, 10)
sol = sympy.solve(e_2 - e_1, x, 0.024)
num = float(sol[0])
print(num)

您的代码中没有定义 f_x 和 g_x 的值。
f_x 和 g_x 是什么？如果您的意思是 e_1 和 e_2，则该方程根本无法解析解。
使用 \"nsolve\" 而不是 \"solve\" 检索数值解。
diameter 未定义。请在发布之前在新的 python 进程中实际测试代码。
为什么要使用 ** 运算符？读起来很糟糕，执行起来效率很低。只需除以价值。

标签： python math sympy equation-solving

【解决方案1】：

通常，nsolve 是用于数值求解方程（或方程组）的 SymPy 工具。但是，我无法使用：它不断引发错误。问题是您的函数是在一个非常小的区域上定义的，并且零非常接近边界：

所以，在这种情况下，我们可以试试numerical root finding techniques。同样，某些工具可能会失败，但经过几次尝试后，我发现 bisect 工作正常：

from sympy import *
from scipy.optimize import root, brentq, bisect
x = symbols('x')

# you didn't provide the diameter, so I've computed it based on your expected solution
d = 1.56843878182221
e_1 = x**-0.5
e_2 = -2 * log(0.00013 * 7-1*d-1+2.51350000**-1*x**-0.5, 10)

# convert the symbolic expression to a numerical function
ff = lambdify(x, e_1 - e_2)
root, output = bisect(ff, 0.023, 0.025, full_output=True)
print(root)
# 0.024000000001862646
print(output)
#      converged: True
#           flag: 'converged'
# function_calls: 32
#     iterations: 30
#           root: 0.024000000001862646

【讨论】：

【解决方案2】：

fixed point method 非常适合用于此类情况。或者，至少将方程转换为兼容形式的原理可以通过提供一种行为不那么恶劣的函数形式来使标准求解器受益。

你有一个定义不明确的方程，格式为y - g(y)，其中y = 1/sqrt(x)。所以让我们得到g 的倒数（称之为G），这样我们就可以解决G(y) - G(g(y)) = G(y) - y 了。

>>> g = e_2.subs(1/x**.5, y)
>>> d = Dummy()
>>> G = solve(g - d, y)[0].subs(d, y)
>>> nsolve(G - y, 6)
6.45497224367903
>>> solve(1/x**.5 - _, dict=True)
{x: 0.024}

将方程 f(x) 重新排列为 x - g(x) 形式的过程可能会使用 SymPy 中的内置方法，但通过将每个 x 替换为虚拟变量，求解它，然后再次用x 替换虚拟符号。不同的g 将更有利于找到不同的根，如下例所示，紫色虚线有利于在 1 附近找到根，而实心蓝色则更好地靠近较小的根。

这是“定点形式”功能的一种可能性：

def fixedpoint_Eqs(eq, x=None):
    """rearrange to give eq in form x = g(x)"""
    f = eq.free_symbols
    fp = []
    if x is None:
        assert len(f) == 1, 'must specify x in this case'
        x = list(f)[0]
    Xeq = eq.replace(lambda _:_ == x, lambda _:Dummy())
    X = Xeq.free_symbols - f
    reps = {xi: x for xi in X}
    for xi in X:
        try:
            g = solve(Xeq, xi, dict=True)
            if len(g) != 1:
                raise NotImplementedError
            fp.append(Eq(x, g[0][xi].xreplace(reps)))
        except NotImplementedError:
            pass
    return fp

>>> fixedpoint_Eqs(x+exp(x)+1/x-5)
Eq(x, -1/(x + exp(x) - 5))
Eq(x, -exp(x) + 5 - 1/x)
Eq(x, log(-x + 5 - 1/x))

【讨论】：