Całkowanie metodą Monte Carlo

0

Cześć,
Czy mógłby ktoś sprawdzić mój program? Docelowo ma liczyć całki.

import numpy as np
s=int(input("szczaly"))
def f(x):
    return x*x+2*x
xk=float(input("xk(koniec calkowania)= "))
xp=float(input("xp(poczatek calkowania)= "))
traf=0
dx=np.abs(xk-xp)
for i in range(s):
    x=xp+np.random.random()*(abs(xp-xk))
    traf+=f(xp+x)
wynik=dx*traf/s
print(wynik)

Z góry dzięki

2

Po modyfikacji jednej linii działa poprawnie dla xp <= xk.
Na dole obliczenia całki rzeczywistej i numerycznie dla losowych przedziałów <xp, xk>, w ostatniej kolumnie błąd względny w procentach.


import numpy as np

def go():
    s=int(input("szczaly"))
    xk=float(input("xk(koniec calkowania)= "))
    xp=float(input("xp(poczatek calkowania)= "))
    print(mc(s, xk, xp))

def f(x):
    return x*x+2*x

def antiderivative(x):
    return x ** 3 / 3. + x ** 2

def mc(s, xk, xp):
    traf=0
    dx=np.abs(xk-xp)    
    for i in range(s):
        x=xp + np.random.random()*(abs(xp-xk))
        #traf+=f(xp+x)
        traf+=f(x)
    wynik=dx*traf/s    
    return wynik

def do_some_testing():
    for i in xrange(10):
        xp = np.random.randint(-100, 100)
        xk = np.random.randint(-100, 100)    
        
        xp, xk = sorted([xp, xk])
        
        numerical = mc(100000, xk, xp)
        real = antiderivative(xk) - antiderivative(xp)
        
        print "{:<4d} xp={:=3d} xk={:=3d} real={:>10.2f} numerical={:>10.2f} error={:>11.8f}". \
            format(i, xp, xk, real, numerical, 100. * (real - numerical) / (real or 1))

do_some_testing()
>>> python mc.py
0    xp=-96 xk= 19 real= 288343.33 numerical= 289232.05 error=-0.30821547
1    xp=-12 xk= 85 real= 212365.33 numerical= 210084.50 error= 1.07401561
2    xp=-25 xk=-10 real=   4350.00 numerical=   4334.93 error= 0.34651047
3    xp=-56 xk=-12 real=  54970.67 numerical=  54896.42 error= 0.13506316
4    xp=-31 xk=  0 real=   8969.33 numerical=   8935.76 error= 0.37426187
5    xp=-27 xk=-17 real=   4483.33 numerical=   4480.90 error= 0.05421452
6    xp= 59 xk= 67 real=  32802.67 numerical=  32802.43 error= 0.00071831
7    xp=-78 xk= 29 real= 161070.67 numerical= 161489.93 error=-0.26029630
8    xp=-50 xk= 23 real=  43751.33 numerical=  43732.98 error= 0.04194458
9    xp= 65 xk= 98 real= 227568.00 numerical= 227536.61 error= 0.01379191
0

Dzięki

1 użytkowników online, w tym zalogowanych: 0, gości: 1