Python Use Python to integrate a function constructed in Mathematica

[member 428835]

Hi PF!

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

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!

 

[willem2]

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

 

[member 428835]

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:

(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)))

 

[willem2]

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

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))

 

[pbuk]

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?

 

[member 428835]

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.