Całkowanie Monte Carlo w Pythonie

Czy kod dla całki podwójnej i potrójnej jest poprawny?
Wyniki wychodzą poprawne ale już przy całkach 2D i 3D

import random
import numpy as np

def fx(x):
    return 4.0/ (1.0+x**2)

def fxy(x,y):
    return x+y

def fxyz(x,y,z):
    return x+y+z

def mc_calka(fx,a,b,N):
    
    hx= float(b-a)/float(N)
    
    arx=np.zeros(N)
    
    for i in range(len(arx)):
        arx[i]=random.uniform(a,b)
    
    integralx=0.0    
    for i in arx:
        integralx+= fx(i)
    
    ans= hx*integralx
    return ans

def mc_calka_double(fxy,a,b,c,d,N):
    
    hx=float((b-a)) / float(N)
    hy=float((d-c)) / float(N)
    
    arx=np.zeros(N)
    ary=np.zeros(N)
    for i in range(len(arx)):
        arx[i]= random.uniform(a,b)
        ary[i]= random.uniform(c,d)
        
    integralxy = 0.0
    for x in arx:
        for y in ary:
            integralxy += fxy(x,y)
    ans = hx*hy * integralxy
    return ans

def calka_triple(fxyz,a,b,c,d,e,f,N):
    hx=float((b-a)) / float(N)
    hy=float((d-c)) / float(N)
    hz=float((f-e)) / float(N)
    
    arx=np.zeros(N)
    ary=np.zeros(N)
    arz=np.zeros(N)
    for i in range(len(arx)):
        arx[i]= random.uniform(a,b)
        ary[i]= random.uniform(c,d)    
        arz[i]= random.uniform(e,f)
        
    integralxyz = 0.0
    for x in arx:
        for y in ary:
            for z in arz:
                integralxyz += fxyz(x,y,z)
                
    ans = hx*hy*hz * integralxyz
    return ans
            
a=0
b=1
c=0
d=1
e=0
f=1
N=100
c1=mc_calka(fx,a,b,N)
print(c1)

c2=mc_calka_double(fxy,a,b,c,d,N)
print(c2)

c3=calka_triple(fxyz,a,b,c,d,e,f,N)
print(c3)

Zakładam, że chodzi o to, że coraz wolniej działa program. To dlatego, że mając N=100, przy liczeniu całki trzeciego stopnia robisz 100^3 iteracji co już daje milion operacji.

Możesz to uprościć jeżeli przyjmiesz, że zawsze losujesz tylko i wyłącznie 100 punktów. Wtedy hx, hy, hz muszą być mniejsze. Chyba w mianowniku będzie pierwiastek k-tego stopnia z N.

Natomiast możesz też liczyć całkę w inny sposób. Załóżmy, że masz funkcję f(x) = x. Wtedy licząc ją w granicy [0,1] możesz znaleźć jej maksimum w punkcie 1. Liczysz pole (1-0) * 1 = 1. Teraz losujesz N punktów k+1 wymiarowych i sprawdzasz czy x_4 < f(x_1, x_2, x_3). Dla naszej funkcji liniowej będzie prościej, tj. x_2 < f(x_1). Jeżeli x_2 jest mniejsze to zwiększasz licznik bo wylosowałeś punkt pod wykresem funkcji. Na końcu (licznik / N)*pole_kwadratu to będzie Twoja całka. W ten sposób będziesz miał znacznie mniej operacji mnożenia i dzielenia więc nie będziesz się musiał tak obawiać błędów zaokrągleń.

@twoj_stary_pijany:

Czyli chyba coś takiego będzie.
Zdefiniowałem 2 funkcje pomocnicze do sprawdzania warunku.
Czy on będą dobre?

import random
import numpy as np

def fxy(x,y):
    return x+y

def fxyz(x,y,z):
    return x+y+z
    

def gxy(x,y):
    return  (1 if (0 <= x <= 1 and 0 <= y <= 1) else -1) 

def gxyz(x,y,z):
    return 1 if (0 <= x <= 1 and 0 <= y <= 1 and 0 <= z <= 1) else -1

def mc_2(fxy,a,b,c,d,N):
    
    tmp=0.0
    licznik=0.0
    arx=np.random.uniform(a,b,N)
    ary=np.random.uniform(c,d,N)
    
    for i in range(len(arx)):
        for j in range(len(ary)):
            if (gxy(arx[i], ary[j]) >=0):
                licznik+=1
                tmp+= fxy(arx[i],ary[j])
    tmp=tmp/ float(licznik)            
    area= licznik / float(N**2) * (b-a)*(d-c)
    return tmp*area

def mc_3(fxyz,a,b,c,d,e,f,N):
    
    tmp=0.0
    licznik=0.0
    arx=np.random.uniform(a,b,N)
    ary=np.random.uniform(c,d,N)
    arz=np.random.uniform(c,d,N)
    
    for i in range(len(arx)):
        for j in range(len(ary)):
            for k in range(len(arz)):
                if (gxyz(arx[i], ary[j], arz[k]) >=0):
                    licznik+=1
                    tmp+= fxyz(arx[i],ary[j], arz[k])
    tmp=tmp/ float(licznik)            
    volume= licznik / float(N**3) * (b-a)*(d-c)*(f-e)
    return tmp*volume
            
a=0
b=1
c=0
d=1
e=0
f=1
N=100
print('start')
c1=mc_2(fxy,a,b,c,d,N)
print(c1)

c2=mc_3(fxyz,a,b,c,d,e,f,N)
print(c2)