在 python 中使用梯形数值积分

Question

我需要在 python 中使用梯形规则集成此功能：

theta = .518/r^2 * dr/(sqrt(2*1.158 + 2/r - .518^2/2r^2))

我已经编写了我的代码，绘制时我应该会看到一个椭圆体结构。 theta 应该运行从 0 到 2pi 并且 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()

我想我犯了一些非常基本的错误，因为我没有从中得到任何情节。我认为梯形例程是正确的，因为它适用于其他功能。非常感谢您的帮助

Answer 1

如果打印数值和 theta 的长度，您会发现它们是空的 lists/arrays。

尝试以下操作：

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()

Answer 2

这里有几个问题。

第一个是np.arrange中的第三个参数不是要生成的值的数量而是步骤。这意味着 theta 将只有一个值，而 n 因此 intg 将有 10 个而不是 100 个值。

假设这是你的意图（100 个值）你可以这样做

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)

然后你绘制 numerical 这基本上是一个数字，你可能想要绘制的是整数值 intg - 为此你还需要转换 intg 从列表到 np.array:

intg = np.array(intg)

通过这些更改，程序可以按预期运行，

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

Answer 3

或者，您可以使用 quadpy（我的一个项目）。

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)

给出0.96194633532。然而，梯形公式主要是出于历史目的而实施的。更好且同样简单的规则是 quadpy.line_segment.Midpoint。更好的方法当然是自适应正交

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)

它给出了更准确的 0.961715309492，甚至 tanh-sinh 求积

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)

给出 0.9617153094932353183036398697528.

在 python 中使用梯形数值积分

using trapezoidal numerical integration in python

python

numerical-methods

numerical-integration