I need to integrate this function using trapezoidal rule in python:
theta = .518/r^2 * dr/(sqrt(2*1.158 + 2/r - .518^2/2r^2))
I have written my code and I should be seeing an ellipsoidal structure when plotted. theta should run from 0 to 2pi and r_min = .16 & r_max = .702
import numpy as np
import matplotlib.pyplot as plt
def trapezoidal(f, a, b, n):
h = float(b-a)/n
result = 0.5*f(a) + 0.5*f(b)
for i in range(1, n):
result += f(a + i*h)
result *= h
return result
intg =[]
v = lambda r: (0.5108/(r**2))* (1./np.sqrt(2*1.158+(2/r)-.5108**2/(2*r**2)))
n = np.arange(1,1000,100)
theta = np.arange(0,2*np.pi,100)
for j in n:
numerical = trapezoidal(v, .16,.702 , j)
intg.append(numerical)
plt.plot(numerical,theta)
plt.show()
I am doing some very elementary mistake I guess, because I am getting no plot out of it. I think the trapezoidal routine is correct, because it worked for other functions. your help is very appreciated
Alternatively, you could use quadpy (a project of mine).
import numpy as np
import quadpy
val = quadpy.line_segment.integrate_split(
lambda r: 0.5108/r**2 / np.sqrt(2*1.158 + 2/r - 0.5108**2/(2*r**2)),
0.15, 0.702, 100,
quadpy.line_segment.Trapezoidal()
)
print(val)
gives 0.96194633532. The trapezoidal formula is mostly implemented for historical purposes, however. A better and equally simple rule is quadpy.line_segment.Midpoint. An even better approach is certainly adaptive quadrature
val, error_estimate = quadpy.line_segment.integrate_adaptive(
lambda r: 0.5108/r**2 / np.sqrt(2*1.158 + 2/r - 0.5108**2/(2*r**2)),
[0.15, 0.702],
1.0e-10
)
print(val)
which gives the more accurate 0.961715309492, or even tanh-sinh quadrature
val, error_estimate = quadpy.line_segment.tanh_sinh(
lambda r: 0.5108/r**2 / np.sqrt(2*1.158 + 2/r - 0.5108**2/(2*r**2)),
0.15, 0.702,
1.0e-30
)
print(val)
which gives 0.9617153094932353183036398697528.