# MaxDeriv4 findet Maximum der 4. Ableitung von f
def MaxDeriv4(f,a,b):
q = derivative(f(x),x,4)
v = q.find_maximum_on_interval(a,b)
return v[0]
# Restglied
def Rf(h,a,b,f4max):
return (b - a)/2880*h*h*h*h*f4max
|
|
# Berechnet Anzahl der nötigen Teilstücke für einen Fehler < threshold
def getN(f,a,b,threshold):
f4max = MaxDeriv4(f,a,b)
N = 1000
for i in range(1, 1000, 10):
h = (b-a)/i
if (abs(Rf(h,a,b,f4max)) < threshold):
N = i
return N
|
|
# Fassregel
def fass(f,n,a,b):
h = (b-a)/n
sum1 = 0
for j in range(1,n):
sum1 = sum1 + f(a+j*h)
sum2 = 0
for j in range(1,n+1):
sum2 = sum2 + f((a+j*h + a+(j-1)*h)/2)
return h/3*( (f(a)+f(b))/2 + sum1 + sum2 + sum2)
|
|
# Substituieren, um nur von 0 bis 1 numerisch integrieren zu müssen
f(x) = 1/(1+x^2)
g(y) = f(x).subs(x=y/(1-y)) * derivative(y/(1-y))
g2(y) = g(y).full_simplify()
print g2
plot(g2,-0.2,1.2)
y |--> 1/(2*y^2 - 2*y + 1) y |--> 1/(2*y^2 - 2*y + 1) |
N = getN(g2,0,1,10e-11)
print N
result = fass(g2, N, 0, 1)
print result.n()
231 1.57079632679490 231 1.57079632679490 |
abs(pi/2 - result).n(digits=20)
1.3058537541230097467e-16 1.3058537541230097467e-16 |
|
|
|
|