在短时间内求解包含积分表达式的方程
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 几乎没有时间就做到了,所以可能存在根本性的错误。
我正在尝试使用 SymPy 求解以下包含积分的方程:
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 几乎没有时间就做到了,所以可能存在根本性的错误。