Solve an equation containing Integral expression in short time

Question

I am trying to solve the following equation containing integrals using SymPy: I have tried to calculate just the integral part using the code below, but it takes a long time to generate expressions in 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()

It takes enormous amount of time to solve the equation. How can i solve the entire equation in short time using SymPy or any other package/library (scipy?)

Posting on behalf of my friend.

Answer 1

SymPy seems to be having a hard time doing that integral. I waited for about 2 minutes on my machine and nothing popped up. Maybe it's not solvable analytically.

So I went for a numerical approach using 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)

Which produces the following output for me:

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

The root can also be seen visually using a plot on MatPlotLib or Desmos:

I assume the time taken is reasonable since it had to evaluate 14 pretty difficult integrals. However, Desmos did it in almost no time so maybe there is something fundamentally wrong.