Use Python to integrate a function constructed in Mathematica

• member 428835
I haven't done much in python. Can I just replace the 180 with pi?Yes, you can replace the 180 with pi. In python, math.pi represents the value of pi. So you can use math.pi instead of manually typing 180.
member 428835
Hi PF!

I followed someone's help on here and have the following code in python that performs monte carlo integration

Python:
from math import *
from random import *

def integrate(alpha):
#  MONTE CARLO INTEGRATION OVER NON-RECTANGULAR DOMAINS
def f(pt):
# RETURN INTEGRAND AS FUNCTION
x = pt[0]
y = pt[1]
return x*y**2

def testRegion(pt):
return (pt[0]/2 > pt[1])

def genpoint():
# GENERATE COORDINATES IN A SQUARE
x = (1 - math.sin(alpha))*random() + math.sin(alpha)
y = (1 - math.sin(alpha))*random() + math.sin(alpha)
return (x,y)

# INITIALIZE
Sum = 0.0
Area = (1 - math.sin(alpha))**2
N = 0

samp_pts = 10000
int_dist = []
iterations = 1000
# INTEGRATION
for _ in range(iterations):
for i in range(samp_pts):
N+=1
pt = genpoint()
if testRegion(pt):
Sum += f(pt)

sol = Sum*Area/N
int_dist.append(sol)

return int_dist

However, in line 10 we see the integrand is typed by a user. Is there a way to import a Mathematica function into python to automate the procedure? Thanks so much for your help!

joshmccraney said:
However, in line 10 we see the integrand is typed by a user. Is there a way to import a Mathematica function into python to automate the procedure? Thanks so much for your help!
You can use Export["filename.txt", expression, "PythonExpression"] in mathematica to write to a textfile, and then read the file in python and use eval() in python to compute the results of the expressions.

https://reference.wolfram.com/language/ref/Export.html?view=all

Thanks for this. However, I'm getting the error "could not be encoded as a Python literal expression." Any ideas? The function looks like this:
Code:
(1/((1 + y)^2))1.60617 (1 - y)^2 (4 + y) BesselJ[4,
304.689 Sqrt[1 - x^2]] BesselJ[4, 304.689 Sqrt[1 - y^2]] Cosh[
5.31755 (1 + (1 - x) Csc[\[Pi]/180])] Cosh[
5.31755 (1 + (1 - y) Csc[\[Pi]/180])] (-(((4 - x) (1 + x)^2)/(
120 (1 - x)^2)) + ((1 - x)^2 (4 +
x) (((4 - Cos[\[Pi]/180]) (1 + Cos[\[Pi]/180]))/(
60 (1 - Cos[\[Pi]/180])^2) - (1 + Cos[\[Pi]/180])^2/(
120 (1 - Cos[\[Pi]/180])^2) + ((4 - Cos[\[Pi]/180]) (1 +
Cos[\[Pi]/180])^2)/(
60 (1 - Cos[\[Pi]/180])^3) - ((4 - Cos[\[Pi]/180]) (1 +
Cos[\[Pi]/180])^2 Cot[\[Pi]/180] Csc[\[Pi]/180])/(
120 (1 - Cos[\[Pi]/180])^2)))/((1 +
x)^2 ((1 - Cos[\[Pi]/180])^2/(1 + Cos[\[Pi]/180])^2 - (
2 (1 - Cos[\[Pi]/180])^2 (4 + Cos[\[Pi]/180]))/(1 +
Cos[\[Pi]/180])^3 - (
2 (1 - Cos[\[Pi]/180]) (4 + Cos[\[Pi]/180]))/(1 +
Cos[\[Pi]/180])^2 - ((1 - Cos[\[Pi]/180])^2 (4 +
Cos[\[Pi]/180]) Cot[\[Pi]/180] Csc[\[Pi]/180])/(1 +
Cos[\[Pi]/180])^2)))

I'm Sorry, but it seems PythonExpression doesn't do variables.
You can convert your expression with CForm[expression], and you can insert the result in a python program. You must first define all functions that you might need. I've done all the functions that this one requires. I hope I have the right function for BesselJ
Python:
import math, scipy.special

Pi = math.pi

def Sqrt(x):
return math.sqrt(x)

def Power(a,b):
return a**b

def BesselJ(z,v):
return scipy.special.jv(z,v)

def Cosh(x):
return math.cos(x)

def Cos(x):
return math.cosh(x)

def Csc(x):
return 1/math.sin(x)

def Cot(x):
return 1/math.tan(x)

def test(x,y):
return (1.60617*Power(1-y,2)*(4+y)*BesselJ(4,304.689*Sqrt(1-Power(x,2)))*BesselJ(4,304.689*Sqrt(1-Power(y,2)))*Cosh(5.31755*(1+(1-x)*Csc(Pi/180.)))*Cosh(5.31755*(1+(1-y)*Csc(Pi/180.)))*(-0.008333333333333333*((4-x)*Power(1+x,2))/Power(1-x,2)+(Power(1-x,2)*(4+x)*(((4-Cos(Pi/180.))*(1+Cos(Pi/180.)))/(60.*Power(1-Cos(Pi/180.),2))-Power(1+Cos(Pi/180.),2)/(120.*Power(1-Cos(Pi/180.),2))+((4-Cos(Pi/180.))*Power(1+Cos(Pi/180.),2))/(60.*Power(1-Cos(Pi/180.),3))-((4-Cos(Pi/180.))*Power(1+Cos(Pi/180.),2)*Cot(Pi/180.)*Csc(Pi/180.))/(120.*Power(1-Cos(Pi/180.),2))))/(Power(1+x,2)*(Power(1-Cos(Pi/180.),2)/Power(1+Cos(Pi/180.),2)-(2*Power(1-Cos(Pi/180.),2)*(4+Cos(Pi/180.)))/Power(1+Cos(Pi/180.),3)-(2*(1-Cos(Pi/180.))*(4+Cos(Pi/180.)))/Power(1+Cos(Pi/180.),2)-(Power(1-Cos(Pi/180.),2)*(4+Cos(Pi/180.))*Cot(Pi/180.)*Csc(Pi/180.))/Power(1+Cos(Pi/180.),2)))))/Power(1+y,2)

print(test (0.2, 0.3))

member 428835
Python:
def test(x,y):
return (1.60617*Power(1-y,2)*(4+y)*BesselJ(4,304.689*Sqrt(1-Power(x,2)))*BesselJ(4,304.689*Sqrt(1-Power(y,2)))*Cosh(5.31755*(1+(1-x)*Csc(Pi/180.)))*Cosh(5.31755*(1+(1-y)*Csc(Pi/180.)))*(-0.008333333333333333*((4-x)*Power(1+x,2))/Power(1-x,2)+(Power(1-x,2)*(4+x)*(((4-Cos(Pi/180.))*(1+Cos(Pi/180.)))/(60.*Power(1-Cos(Pi/180.),2))-Power(1+Cos(Pi/180.),2)/(120.*Power(1-Cos(Pi/180.),2))+((4-Cos(Pi/180.))*Power(1+Cos(Pi/180.),2))/(60.*Power(1-Cos(Pi/180.),3))-((4-Cos(Pi/180.))*Power(1+Cos(Pi/180.),2)*Cot(Pi/180.)*Csc(Pi/180.))/(120.*Power(1-Cos(Pi/180.),2))))/(Power(1+x,2)*(Power(1-Cos(Pi/180.),2)/Power(1+Cos(Pi/180.),2)-(2*Power(1-Cos(Pi/180.),2)*(4+Cos(Pi/180.)))/Power(1+Cos(Pi/180.),3)-(2*(1-Cos(Pi/180.))*(4+Cos(Pi/180.)))/Power(1+Cos(Pi/180.),2)-(Power(1-Cos(Pi/180.),2)*(4+Cos(Pi/180.))*Cot(Pi/180.)*Csc(Pi/180.))/Power(1+Cos(Pi/180.),2)))))/Power(1+y,2)

What is the point of repeatedly evaluating all those constants like ## 1 - \cos \frac{\pi}{180} ##? And you do realize that both Mathematica and Python work in radians and so 180 is completely meaningless in this context?

What are you actually trying to achieve?

Last edited by a moderator:
pbuk said:
Python:
def test(x,y):
return (1.60617*Power(1-y,2)*(4+y)*BesselJ(4,304.689*Sqrt(1-Power(x,2)))*BesselJ(4,304.689*Sqrt(1-Power(y,2)))*Cosh(5.31755*(1+(1-x)*Csc(Pi/180.)))*Cosh(5.31755*(1+(1-y)*Csc(Pi/180.)))*(-0.008333333333333333*((4-x)*Power(1+x,2))/Power(1-x,2)+(Power(1-x,2)*(4+x)*(((4-Cos(Pi/180.))*(1+Cos(Pi/180.)))/(60.*Power(1-Cos(Pi/180.),2))-Power(1+Cos(Pi/180.),2)/(120.*Power(1-Cos(Pi/180.),2))+((4-Cos(Pi/180.))*Power(1+Cos(Pi/180.),2))/(60.*Power(1-Cos(Pi/180.),3))-((4-Cos(Pi/180.))*Power(1+Cos(Pi/180.),2)*Cot(Pi/180.)*Csc(Pi/180.))/(120.*Power(1-Cos(Pi/180.),2))))/(Power(1+x,2)*(Power(1-Cos(Pi/180.),2)/Power(1+Cos(Pi/180.),2)-(2*Power(1-Cos(Pi/180.),2)*(4+Cos(Pi/180.)))/Power(1+Cos(Pi/180.),3)-(2*(1-Cos(Pi/180.))*(4+Cos(Pi/180.)))/Power(1+Cos(Pi/180.),2)-(Power(1-Cos(Pi/180.),2)*(4+Cos(Pi/180.))*Cot(Pi/180.)*Csc(Pi/180.))/Power(1+Cos(Pi/180.),2)))))/Power(1+y,2)

What is the point of repeatedly evaluating all those constants like ## 1 - \cos \frac{\pi}{180} ##? And you do realize that both Mathematica and Python work in radians and so 180 is completely meaningless in this context?

What are you actually trying to achieve?
This is the product of two basis functions ##\phi## and a Green's function. The basis functions satisfy ##\nabla^2 \phi = 0## in the domain and and ##\phi_n = 0## on the boundary. Then it's evaluated along an arc. I do this all symbolically for the geometry of interest, so I'm not actually typing this stuff out. However, to get something PF can reiterate, I have to evaluate it.

At the end of the day I'm integrating this function over a triangular domain. Sometimes that integral has issues, so someone recommended I do Monte Carlo integration. The native MCI to Mathematica takes FOREVER and is sometimes inaccurate, so decided to make my own in python, per the guidance of another user. Issue here is importing this function into python. But, turns out it's not too bad.

And the 180 is not meaningless: the angle I specified was something like ##89 \pi/180##, so something like that. Degrees are easier for me to track than radians.

1. How can I use Python to integrate a function constructed in Mathematica?

Python has a built-in function called scipy.integrate.quad that can be used to integrate any function, including those constructed in Mathematica. This function takes in the name of the function, the lower and upper limits of integration, and any additional parameters. It returns both the value of the integral and an estimate of the error.

2. Can I pass a Mathematica function directly to Python?

Yes, you can pass a Mathematica function to Python by first exporting it as a Python module using the Export function in Mathematica. Then, you can import the module in Python and use the function just like any other Python function.

3. Is there a performance difference between integrating a function in Mathematica vs. Python?

In general, integrating a function in Mathematica is faster than using Python's scipy.integrate.quad function. However, the difference in performance may vary depending on the complexity of the function and the specific parameters used.

4. Can I use Python libraries within a Mathematica function?

Yes, you can use Python libraries within a Mathematica function by using the ExternalEvaluate function. This allows you to call Python code and pass data between Mathematica and Python.

5. Are there any limitations to integrating a Mathematica function in Python?

There are a few limitations to keep in mind when integrating a Mathematica function in Python. One limitation is that the function must be numerical, meaning it cannot contain symbolic or complex expressions. Additionally, the function must be defined using standard mathematical operators and functions that are supported by both Mathematica and Python.

• Programming and Computer Science
Replies
8
Views
845
• Programming and Computer Science
Replies
15
Views
1K
• Programming and Computer Science
Replies
55
Views
4K
• Programming and Computer Science
Replies
1
Views
1K
• Programming and Computer Science
Replies
1
Views
2K
• Programming and Computer Science
Replies
1
Views
1K
• Programming and Computer Science
Replies
1
Views
707
• Programming and Computer Science
Replies
18
Views
1K
• Programming and Computer Science
Replies
16
Views
1K
• Programming and Computer Science
Replies
7
Views
3K