我试图通过复合梯形规则与通过精确公式来估计计算高斯积分的百分比误差。为此，我在 Python 中生成了一个平均值为 0、标准偏差和振幅范围内个点的网格上, 并用于计算一系列不同的近似积分$\sigma = 1$$A = 1$$N$$[-2\sigma,2\sigma]$scipy.integrate.trapz$A$$N$

对于精确积分，我使用了公式计算两个均值标准差之间的精确积分，因为，或$I = \int_{-x}^{x}e^{\frac{-t^2}{2\sigma^2}} dt = \sqrt{2\pi}\sigma\text{erf}(\frac{x}{\sqrt{2}\sigma})$$\pm$$I = \int_{-2\sigma}^{2\sigma}e^{\frac{-t^2}{2\sigma^2}} dt = \sqrt{2\pi}\sigma\text{erf}(\sqrt{2})$

$I = \sqrt{2\pi}\text{erf}(\sqrt{2})$

在这种的情况下。那么 % 错误是。下面是与$\sigma = 1$$E = \frac{|I - A|}{I}\cdot 100$$\frac{1}{E}$$N$

并作为参考 verus的图。$E$$N$

现在我希望但这在这里看不到。谁能帮我弄清楚为什么？$E \propto \frac{1}{N^2}$

我在下面包含了我的代码。先感谢您！

from scipy.integrate import trapz
from scipy.special import erf
import matplotlib.pyplot as plt
import numpy as np

errs = np.array([]) # Errors in integral estimates
Ns = np.array([]) # Sample sizes tried

def get_percent_err(N, errs, Ns):
    A = 1 # Amplitude of Gaussian
    std = 1 # Standard deviation of Gaussian
    r = 2*std # Extent of domain
    dt = 2*r/N # Step
    t = np.linspace(-r,r,N) # Domain
    gauss = A*np.exp(-((t)**2)/(2*(std**2))) # Gaussian function
    
    actual_area = np.sqrt(2*np.pi)*A*std*erf(np.sqrt(2)) # Exact gaussian integral from -2*std to +2*std
    est_area = trapz(gauss, dx=dt) # Estimated integral via trapezium rule
    percent_err = (abs(actual_area - est_area)/actual_area)*100 # Percentage error
    
    errs = np.append(errs, percent_err)
    Ns = np.append(Ns, N)
    
    return errs, Ns

# Sweep through a range of sample sizes and obtain errors for each
min_samples = 100
max_samples = 100100
sample_step = 100
samples = np.arange(min_samples,max_samples,sample_step)

for N in samples:
    errs, Ns = get_percent_err(N, errs, Ns)

# Plot reciprocal of error in integral against sample size
plt.figure()
plt.plot(1/Ns,errs)
plt.xlabel("1 / Samples (1 / N)")
plt.ylabel("% Error in Integral (E)")
plt.title("Plot of % Error (E) in Integral versus 1 / Samples (1 / N)")

更新：

由于我不知道的原因（尽管我相信它与数值精度有关），此代码改编自此链接，给出了预期的结果。

代码：

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erf
from scipy.integrate import trapz as trapz_scipy
from pandas import DataFrame

def trapz(f,a,b,N=50):
    x = np.linspace(a,b,N+1) # N+1 points make N subintervals
    y = f(x)
    dx = (b - a)/N
    T = trapz_scipy(y, dx=dx)
    return T

def get_err(N, errs):
    A = 1
    std = 1
    f = lambda x : A*np.exp(-x**2/(2*std**2))
    a = -2*std; b = 2*std;
    
    I = np.sqrt(2*np.pi)*A*std*erf(np.sqrt(2))
    T = trapz(f,a,b,N)
    return np.append(errs, 100*np.abs(I - T)/T)

min_samples = 10
max_samples = 1000
sample_step = 10
sample_sizes = np.arange(min_samples,max_samples+sample_step,sample_step)

errs = np.array([], dtype=np.float64)

for s in sample_sizes:
    errs = get_err(s, errs)
    
err_vs_sample_size = DataFrame({'% Error' : errs, 'Samples' : sample_sizes})

plt.figure()
plt.plot(1/sample_sizes**2,errs)
plt.xlabel("1 / Samples (1 / N)")
plt.ylabel("% Error in Integral (E)")
plt.title("Plot of % Error (E) in Integral versus 1 / Samples (1 / N)")

错误图：