数值积分问题

Question

我目前被一个隐藏在某个地方的非常简单的错误所困，希望有人能对此有所了解。

我正在尝试将一组数据与常规方法进行数值整合，并获得我不太期望的结果。

然后我回到第一原则并生成了以下代码。由于某种原因，我无法生成 y3 或 y4 来匹配 y2。

import numpy as np  
import matplotlib.pyplot as plt
import math

f1 = 'x**2 - 1' # Starting function
start, stop, step = -2, 2, 21

x = np.linspace(start, stop, step)
xdash = []
xDoubleDash = []

y1 = eval(f1)
y2 = [] # integrated function
y3 = [] # box rule
y4 = [] # trapezium rule
y5 = [] # numerical differentialtion of integrated function

for i in range(0,len(x)-1):
    xdashElement = (x[i] + x[i+1])/2
    y2.append((math.pow(xdashElement,3)/3) - xdashElement) # Integrated Function: (x**3)/3 - x
    xdash.append(xdashElement)
    y3.append((y1[i]+y1[i+1])/2 * abs(((x[i] - x[i+1])))) # box rule
    y4.append(((y1[i]+y1[i+1]) * abs(x[i] - x[i+1]))/2) # trapezium rule

for i in range(0,len(y2)-1):
    xDoubleDashElement = (xdash[i] + xdash[i+1])/2
    xDoubleDash.append(xDoubleDashElement)
    y5.append((y2[i+1]-y2[i])/(xdash[i+1] - xdash[i])) 


plt.plot(x, y1, 'b-')
plt.plot(xdash, y2, 'r-')
plt.plot(xdash, y3, 'g-')
plt.plot(xdash, y4, 'm-')
plt.plot(xDoubleDash, y5, 'c-')
plt.grid()
plt.show()

提前致谢。

编辑：我添加了 y5，y2 的数值微分作为健全性检查。

Answer 1

您的y2 是不定积分，而您的y3 和y4 是一系列区间上的定积分。我不会指望他们彼此匹配。如果你想用反导数检查你的数值积分，你的y3 和y4 应该匹配这个y2：

ya = (math.pow(x[i],3)/3) - x[i]
yb = (math.pow(x[i+1],3)/3) - x[i+1]
y2.append(yb-ya) # Integrated Function: (x**3)/3 - x

编辑以响应integral上的cmets：

我期望 y3 和 y4 的值与 y2 的值​​非常相似，并且如果增量 x[i+1] - x[i] 根据积分的理论定义变得无限小，理论上会匹配

p>

基于this definition of integral，

其中F 是f 的反导数或不定积分。

对于您的 y3 或 y4，x=(x[i] + x[i+1])/2 您有 y：

对于你的y2，x=(x[i] + x[i+1])/2 你有 y：

当“x[i+1] - x[i] 变得无限小”时，你会有

换句话说，当“x[i+1] - x[i] 变得无限小”时，您的y3 或y4 几乎是y=0 处的一条水平线。为什么你的y2 应该匹配y=0 处的一行？

如果你真的想看到你的y3 和y4 匹配y2 的形状，你需要一个y0 或“arbitrary constant of integration”来开始并制作y3 和y4像这样：

y0 = -.7
y3 = [y0] # box rule
y4 = [y0] # trapezium rule
for i in range(0,len(x)-1):
    xdashElement = (x[i] + x[i+1])/2
    y3.append((y1[i]+y1[i+1])/2 * abs(((x[i] - x[i+1]))) + y3[len(y3) - 1]) # box rule
    y4.append(((y1[i]+y1[i+1]) * abs(x[i] - x[i+1]))/2 + y4[len(y4) - 1]) # trapezium rule
y3.pop(0)
y4.pop(0)

然后随着step 的增加或“增量 x[i+1] - x[i] 变得无限小”，这些y3 和y4 将匹配y2。