在短时间内求解包含积分表达式的方程

Solve an equation containing Integral expression in short time

我正在尝试使用 SymPy 求解以下包含积分的方程: 我尝试使用下面的代码只计算整数部分,但是在 r

中生成表达式需要很长时间
from sympy import *
mean,std =0,1
Q=250
#defining Cumulative distribution function
def cdf(mean,std):
  t,x = symbols('t,x')
  cdf_eqn = (1/(std*sqrt(2*pi)))*exp(-(((t-mean)**2)/(2*std**2)))
  cdf = Integral(cdf_eqn, (t,-oo,x)).doit()
  return cdf
#defining Probability density function
def pdf(mean,std):
  x = symbols('x')
  pdf = (1/(std*sqrt(2*pi)))*exp(-((( (x - mean)**2)/(2*std**2)))).doit()
  return pdf
#multiplying cdf and pdf
r,x = symbols('r,x')
equation = cdf(mean=0,std=1).subs(x,x)*pdf(mean=0,std=1).subs(x,(r + Q -x))
#getting interating equation over the limits [0,r]
final_equation = Integral(equation, (x,0,r))
#solving the equation
final_equation.doit()

求解方程需要花费大量时间。我如何使用 SymPy 或任何其他 package/library (scipy?)

在短时间内求解整个方程

代表朋友发帖。

SymPy 似乎很难完成这个积分。我在我的机器上等待了大约 2 分钟,但没有弹出任何内容。也许它无法解析解决。

所以我采用了 SciPy's root finding algorithm 的数值方法。

import sympy as sp
from scipy.optimize import root_scalar
import time

start_time = time.time()

mean, std = 0, 1
Q = 250
p = 5
w = 2

x = sp.symbols("x")
r_symbol = sp.symbols("r")

pdf = (1 / (std * sp.sqrt(2 * sp.pi))) * sp.exp(-(((x - mean) ** 2) / (2 * std ** 2)))
cdf = sp.erf(sp.sqrt(2) * x / 2) / 2 + 1 / 2  # pre-calculated with integrate(pdf, (x, -oo, t))


def f(r: float) -> float:
    result = sp.N(-p + (p + w * cdf.subs(x, Q)) * cdf.subs(x, r) + \
                  w * sp.Integral(cdf * pdf.subs(x, (r + Q - x)), (x, 0, r)))
    return result


r0 = 1  # initial estimate for the root
bracket = (-10, 10)  # the upper and lower bounds of where the root is
solution = root_scalar(f, x0=r0, bracket=bracket)
print(solution)  # info about the convergence
print(solution.root)  # the actual number

end_time = time.time()
print("Time taken:", end_time - start_time)

这会为我生成以下输出:

      converged: True
           flag: 'converged'
 function_calls: 14
     iterations: 13
           root: 0.5659488219328516
0.5659488219328516
Time taken: 26.701611518859863

也可以使用 MatPlotLib 或 Desmos 上的绘图直观地看到根:

我认为花费的时间是合理的,因为它必须评估 14 个非常困难的积分。然而,Desmos 几乎没有时间就做到了,所以可能存在根本性的错误。