0
votes

I have been having some trouble getting the sympy module to evaluate a definite integral. Equation When I try to run the following code the program fails to finish. The problem seems to come from the fact that the integral bounds includes a variable that is in the equation. It keeps running and running and is never able to resolve the equation. I was able to get the equation to resolve by using different bounds but that's not what I'm going for. If you have any advice or tips I'd be glad to hear it. Thanks for your help.

import sympy
import numpy as np
from scipy.integrate import quad
import matplotlib.pyplot as plt

sympy.init_printing()
Q, eta, phi, kappa, lamb, beta, T0, T1, deltaT, T = sympy.symbols('Q eta phi kappa lamb beta T0 T1 deltaT T', real=True)


integrand = sympy.integrate(eta**(lamb*beta-1)*sympy.exp(-eta**2), (eta,0,eta))

T = T0 + deltaT*integrand

fT = sympy.lambdify((T0,T1,eta,lamb,beta), T, 'numpy')
1
whats the output? what bounds are breaking it? - PYA
Your notation is flawed. Do not use the same symbol for integration bound and integration variable. - laolux
Prith, I am trying to make a callable function using sympy from this equation. The bounds that are breaking it is the eta in the integration bounds which the equation calls for. Hannebambel, you mention that I shouldn't be using the same symbol for integration but the equation I am trying to replicate does use the same symbol for integration bound and integration variable. - DanFrankenstein
Yeah, I have seen the equation and it is bad to begin with. Imagine I want to calculate T(2*eta). Where would the 2 appear on the right hand side? Naturally one would expect it appears before every eta on the right hand side. But this would not change the integral at all, so T(eta)=T(2*eta)... - laolux
Ahh, I see what you mean. Well, I'll give that a try. Thanks for the help! - DanFrankenstein

1 Answers

0
votes

Seems like numpy does not have the required functionality for this expression. See my example code below.

import sympy as sym
import numpy as np

eta, eta_prime, lamb, beta, temp = sym.symbols('eta eta_prime lamb beta temp', real=True)
integrand = sym.integrate(eta_prime**(temp)*sym.exp(-eta_prime**2), (eta_prime, 0, eta))
integrand = integrand.replace(temp, lamb*beta-1)
integrand = sym.simplify(integrand)

sym.pprint(integrand)

T, T0, T1, deltaT = sym.symbols('T T0 T1 deltaT', real=True)
T = T0 + deltaT*integrand

sym.pprint(T)

fT = sym.lambdify((T0,T1,eta,lamb,beta),T)

print(fT(1,2,3,4,5))

fT_numpy = sym.lambdify((T0,T1,eta,lamb,beta),T, 'numpy')

print(fT_numpy(1,2,3,4,5))

The last line breaks because numpy does not know about the lowergamma function used by sympy.

Note I had to introduce a temporary real variable temp in the integral which is later replaced by the correct expression lamb*beta-1. Mathematically this replacement is exact, but somehow sympy does not see it and takes forever when given the original expression. I also replaced eta with eta_prime in the integral (not the bounds).