How to solve integrating an equation by python

Question

I am trying to integrate an equation by python. However, I don't understand why the integration doesn't run. The equation I want to integrate is:

Having: , to find having and , with .

I am doing this procedure:

from sympy import *

x=symbols('x')
y=symbols('y')

Gamma = 0.167
sigma_8 = 0.9
M_8 = 6e14

gamma = (0.3*Gamma+0.2)*(2.92+1/3*log(x/M_8))
sigma = 0.9*(x/M_8)**(-gamma/3)
diff(sigma,x)
print(diff(sigma,x))

out:

0.9*(1.66666666666667e-15*x)**(-0.0277888888888889*log(1.66666666666667e-15*x) - 0.243430666666667)*(600000000000000.0*(-4.63148148148148e-17*log(1.66666666666667e-15*x) - 4.05717777777778e-16)/x - 0.0277888888888889*log(1.66666666666667e-15*x)/x)

Then

import math

dnudM = math.sqrt(0.707)*1.686*(1+y)*diff(sigma,x)
print(dnudM)

Out:

0.9*(1.66666666666667e-15*x)**(-0.0277888888888889*log(1.66666666666667e-15*x) - 0.243430666666667)*(1.41764430376593*y + 1.41764430376593)*(600000000000000.0*(-4.63148148148148e-17*log(1.66666666666667e-15*x) - 4.05717777777778e-16)/x - 0.0277888888888889*log(1.66666666666667e-15*x)/x)

Then

n = (1+(1/x**0.6))*exp(-x**2/2)*dnudM
print(n)

And out

0.9*(1.66666666666667e-15*x)**(-0.0277888888888889*log(1.66666666666667e-15*x) - 0.243430666666667)*(x**(-0.6) + 1)*(1.41764430376593*y + 1.41764430376593)*(600000000000000.0*(-4.63148148148148e-17*log(1.66666666666667e-15*x) - 4.05717777777778e-16)/x - 0.0277888888888889*log(1.66666666666667e-15*x)/x)*exp(-x**2/2)

Finally, I arrive to this point that the integration doesn't produce any output.

n_H = integrate(n, x)
print(n_H)

It doesn't show either errors nor output!

Answer 1

The code you have ran to completion for me, but took nearly 30 hours. I tweaked a couple of variable names and added in some timer code to track it.

Code:

import math
from sympy import *
import time
x=symbols('x')
y=symbols('y')

Gamma = 0.167
sigma_8 = 0.9
M_8 = 6e14

gamma_x = (0.3*Gamma+0.2)*(2.92+1/3*log(x/M_8))
sigma = 0.9*(x/M_8)**(-gamma_x/3)
diff_sigma_x = diff(sigma,x)
print(f"{diff_sigma_x=}\n")

dnudM = math.sqrt(0.707)*1.686*(1+y)*diff_sigma_x
print(f"{dnudM=}\n")

n = (1+(1/x**0.6))*exp(-x**2/2)*dnudM
print(f"{n=}\n")

start = time.time()
n_H = integrate(n, x)
print(f"{n_H=}\n")
end = time.time()

print(f"Time to integrate = {end-start} seconds")

Results

n_H=-2.10235303142289*(y + 1)*(1.0*Integral(-4.20001660788705e-11*exp(-x**2/2)/(1.66666666666667e-15**(0.0277888888888889*log(x))*x**0.897831723570746*x**(0.0277888888888889*log(x))), x) + 1.0*Integral(-4.20001660788705e-11*exp(-x**2/2)/(1.66666666666667e-15**(0.0277888888888889*log(x))*x**0.297831723570746*x**(0.0277888888888889*log(x))), x) + 1.0*Integral(1.41662964847297e-12*exp(-x**2/2)*log(x)/(1.66666666666667e-15**(0.0277888888888889*log(x))*x**0.897831723570746*x**(0.0277888888888889*log(x))), x) + 1.0*Integral(1.41662964847297e-12*exp(-x**2/2)*log(x)/(1.66666666666667e-15**(0.0277888888888889*log(x))*x**0.297831723570746*x**(0.0277888888888889*log(x))), x))

Time to integrate = 107187.83417797089 seconds