using trapezoidal numerical integration in python

Question

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

Nico Schlömer · Accepted Answer · 2017-10-09 10:01:10Z

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.

atru · Accepted Answer · 2017-09-26 17:09:55Z

There are a couple of issues here.

First one is that the third argument in np.arrange is not the number of values to be generated but the step. This means that theta will have only one value and that n and thus intg will have 10 instead of 100 values.

Assuming that was your intention (100 values) you can do this

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,10)
theta = np.arange(0,2*np.pi,2*np.pi/100)
#print theta
for j in n:
    numerical = trapezoidal(v, .16,.702 , j)
    intg.append(numerical)

Then you're plotting numerical which is basically a single number and what you probably wanted to plot was the integral value intg - to do so you also need to convert intg from a list into np.array:

intg = np.array(intg)

With these changes the program works as intended,

plt.plot(intg,theta)
plt.show()

GioR · Accepted Answer · 2017-09-26 17:07:29Z

0

If you print the lengths of your numerical and theta, you will see that they are empty lists/arrays.

Try the following:

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, 1001)
theta = np.linspace(0,2.*np.pi,1000)

for j in n:
    numerical = trapezoidal(v, .16,.702 , j)
    intg.append(numerical)


plt.plot(intg,theta)
plt.show()

answered Sep 26, 2017 at 17:07

GioR

6362 gold badges9 silver badges19 bronze badges

Collectives™ on Stack Overflow

using trapezoidal numerical integration in python

3 Answers 3

Comments

Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

3 Answers 3

Comments

Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related