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