函数的np.trapz

Question

我搞不清楚np.trapz。我应该自己编写一个梯形规则，然后将其与np.trapz进行比较。然而，这里有一个陷阱。假设积分是从a到b的。我应该找到a=1b=1，a=0b=2，a=0b=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 itteration 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

您没有指定任何关于误差/公差接受的内容，但在比较数值计算积分的不同方法（..、方程或类似方法）时，通常会比较它们在系统大小、迭代步骤数等增加时的误差比例

梯形法则是一种常见的教科书方法，有大量的自由文献（例如关于Wikipedia）。因此，您可以轻松地验证实现的有效性，因为它的行为已经众所周知。另外，在一个问题上测试你的方法总是一个好主意，这个问题可以通过分析（用笔和纸）来解决。这样做会很快发现代码/实现中的错误

那我们就这么做吧！考虑平凡函数^ {< CD5> }。将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)函数在a和b之间创建一个n线性间隔数数组。a我们现在可以称之为

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中的实现进行比较。在阅读了它的documentation之后，我们看到积分是通过

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_list中n_steps的每个值，使用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。这与维基百科文章提供的错误分析是一致的

现在，我建议您尝试重复这个过程，但是用您的原始函数f(x)替换我的平凡的g(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))