函数的 np.trapz

Question

我无法弄清楚 np.trapz。我应该自己编写一个梯形规则，然后将其与 np.trapz 进行比较但是，有一个问题。假设积分是从a到b。我应该找到 a=1 b=1、a=0 b=2、a=0 b=3... a=0 b=10 的积分并绘制这些值。这是我制作的梯形函数的代码：

## trapezoidal ##
import numpy as np
import matplotlib.pyplot as plt

def f(x):
    o1 = 0.3 ## matter density of universe
    o2 = 0.7 ## density of dark enerrgy
    o3 = 0 ## density of radiation
    c = 3*10**3 ## constant
    return c/((o1*(1+x)**3 + o3*(1+x)**2 + o2)**(1/2))

data = [] ## data array

for b in range (0, 10):## for loop going through z2 as b

    ## definitions 
    a = 0
    N = 100
    h = (b-a)/N
    integral = 0

    integral = 0.5*(f(a) + f(b)) ## endpoints

    for i in range(1,N):
        integral += f(a + (i*h))
    integral = integral*h
    integral = integral/(1+b)
    data.append(integral) ## appending each iteration into data array

 print(data)

plt.plot([1,2,3,4,5,6,7,8,9,10],data)
plt.xlabel('z')
plt.ylabel('DA')

这是我对 np.trapz 所做的尝试，尽管我认为我完全错了。 ## np.trapz ## 将 numpy 导入为 np 将 matplotlib.pyplot 导入为 plt

def f(x):
    o1 = 0.3 ## matter density of universe
    o2 = 0.7 ## density of dark enerrgy
    o3 = 0 ## density of radiation
    c = 3*10**3 ## constant
    return c/((o1*(1+x)**3 + o3*(1+x)**2 + o2)**(1/2))

x = np.arange(1,10,1)
##area = []

for i in range (0,10):
    x = [0,i]
    y = f/(1+i)
    area.append(np.trapz(y,x))

print(area)

plt.plot([1,2,3,4,5,6,7,8,9,10],area)
plt.xlabel('z')
plt.ylabel('DA')

Answer 1

我将尝试展示如何实现梯形规则，并将其与

numpy

的实现进行比较。因此，您有一个函数

def f(x)

，您希望将其从

a

集成到

b

。

您没有指定任何有关误差/容差接受的信息，但在比较不同的数值评估积分方法（..、方程或类似方法）时，通常会比较当系统规模、迭代步数和迭代次数增加时它们的误差规模如何。等等。

梯形规则是一种常见的教科书方法，并且有大量免费文献（例如在维基百科上）。因此，您可以轻松验证实现的有效性，因为其行为已经众所周知。此外，在无法通过分析（通过笔和纸）解决的问题上测试您的方法始终是一个好主意。这样做将很快发现代码/实现中的错误。

所以我们就这么做吧！考虑平凡函数

g(x) = x

。将

g(x)

从

0

积分到

10

很容易完成，并得到

(1/2) * (10^2 - 0^2) = 50

。这是Python中的

g(x)

，

def g(x):
  return x

我们现在实现一个函数，用于评估任何函数的积分

func

，使用维基百科上定义的梯形规则（上面的链接），

def my_trapz(func, a, b, n_steps):
  X = np.linspace(a,b,n_steps)
  integral = (func(a)+func(b))/2.0 + sum([func(x) for x in X])
  return integral * (b-a)/n_steps

这里，

func

是一个Python函数（

def

），只接受一个参数。

a

是积分的下限，

b

是积分的上限，而

n_steps

是要在

[a,b]

范围内计算的 x 值的数量。

n_steps

越大，积分越准确。

numpy.linspace(a,b,n)

函数在

n

和

a

之间创建一个由

b

线性间隔的数字组成的数组。我们现在可以打电话，

a = 0.0
b = 10.0
n_steps = 100

integral = my_trapz(g, a, b, n_steps)
print(integral) # outputs 50.4999999..

您会注意到，为

n_steps

输入的值越高，结果就越接近

50

，正如预期的那样。现在，我们想要将我们的实现与

numpy

中的实现进行比较。阅读其文档后，我们看到积分的计算方式为，

X = np.linspace(a,b,n_steps)

integral = np.trapz([g(x) for x in X], dx=(b-a)/n_steps)
print(integral) # outputs 49.9023437.., slightly better than out implementation

在这里，我们直接在函数调用中使用了列表理解。或者，我们也可以这样做

g_values = []
for x in :
  g_values.append(g(x))
integral = np.trapz(g_values, dx=(b-a)/n_steps)

我们可以通过同意一组不同的

n_steps

来比较它们，并研究这两种方法的执行方式：对于列表

n_steps

中的每个

n_steps_list

值，使用

numpy

和我们自己的计算相应的积分方法。然后将积分结果绘制为

n_steps

的函数。也就是说，

a = 0.0
b = 10.0
n_steps_list = [2**n for n in range(3,10)] # =[8, 16, 32, ..., 1024]

integral_np = []
integral_my = []

for n_steps in n_steps_list:
  X = np.linspace(a,b,n_steps)
  dx = (b-a)/n_steps
  integral_np.append(np.trapz([g(x) for x in X], dx=dx))
  integral_my.append(my_trapz(g, a, b, n_steps))

plt.scatter(n_steps_list, integral_np, color='g', label='np')
plt.scatter(n_steps_list, integral_my, color='r', label='my')
plt.xlabel('number of steps (resolution)')
plt.ylabel('Integral result')
plt.legend()
plt.show()

结果图如下所示，

我们看到两种方法都以指数方式收敛到期望值

50

。这与维基百科文章提供的错误分析一致。

现在，我建议您尝试重复此过程，但将我的琐碎

g(x)

替换为您的原始函数

f(x)

并调查会发生什么。

def f(x):
    o1 = 0.3 ## matter density of universe
    o2 = 0.7 ## density of dark enerrgy
    o3 = 0 ## density of radiation
    c = 3*10**3 ## constant
    return c/((o1*(1+x)**3 + o3*(1+x)**2 + o2)**(1/2))