尝试使用 python 编写 e^(-x^2) 在 0 和 1 之间积分的切比雪夫求积解法

Question

我不知道我做错了什么。我最终得到的答案是 0.84965，而它应该更接近 0.746824

import numpy as np
import math

def f(x):
    return np.e**(-x**2)

def chebyshev(a,b,n):
    h=(b-a)/n
    s=0
    r=(b-a)/2
    for i in range(n+1):
        l=i*h
        if i>n/2:
            l=h*(n/2-(i-n/2))
        val=np.cos(np.arcsin(l/r))*r
        s+=f(val)

    return s*h
print(chebyshev(0,1,1000))

Answer 1

如果我使用这些注释第11页中的方程来进行切比雪夫-高斯近似，并执行以下操作：

from math import pi, cos, sqrt, exp


def chebyshev_node(i, n):
    return cos(0.5 * pi * (2 * i - 1) / n)


def chebyshev_solver(f, a, b, n):
    d = (b - a)
    c = 0.5 * pi * d / n

    s = 0.0
    for i in range(1, n + 1):
        cn = chebyshev_node(i, n)
        v = 0.5 * (cn + 1) * d + a
        s += f(v) * sqrt(1 - cn**2)

    return s * c


def efunc(x):
    return exp(-(x**2))


print(chebyshev_solver(efunc, 0, 1, 1000))

它给出 0.7468244140713791，这似乎符合您预期的解决方案。

更新：请注意，以上内容都可以使用 NumPy 向量化为：

import numpy as np


def chebyshev_nodes(n):
    return np.cos(0.5 * np.pi * (2 * np.arange(1, n + 1) - 1) / n)


def chebyshev_solver(f, a, b, n):
    d = (b - a)
    c = 0.5 * np.pi * d / n

    cn = chebyshev_nodes(n)
    v = 0.5 * (cn + 1) * d + a
    return c * np.sum(f(v) * np.sqrt(1 - cn**2))


def efunc(x):
    return np.exp(-(x**2))

print(chebyshev_solver(efunc, 0, 1, 1000))