我听说过,当一个人试图对表格进行数值积分时
和平滑且表现良好(例如,本身不是高度振荡的、非奇异的等),那么将其重写为将有助于准确性
并首先在数值上执行内积分。我看不出有什么理由可以期望它起作用,但是数值方法的准确性又不是很明显。
当然,我知道实际做到这一点的最佳方法是使用一种针对这样的振荡积分进行优化的方法,但为了好奇,假设我限制自己使用一些求积法则。任何人都可以确认或反驳进行这种转换往往会提高积分的准确性吗?和/或指向我解释它的来源?
我听说过,当一个人试图对表格进行数值积分时
和平滑且表现良好(例如,本身不是高度振荡的、非奇异的等),那么将其重写为将有助于准确性
并首先在数值上执行内积分。我看不出有什么理由可以期望它起作用,但是数值方法的准确性又不是很明显。
当然,我知道实际做到这一点的最佳方法是使用一种针对这样的振荡积分进行优化的方法,但为了好奇,假设我限制自己使用一些求积法则。任何人都可以确认或反驳进行这种转换往往会提高积分的准确性吗?和/或指向我解释它的来源?
我不认为这有什么区别。您必须为积分选择足够高的正交所以它等于贝塞尔函数. 我在下面的示例中选择了 20 阶,但是您始终必须针对您积分的确切函数和区间进行收敛。然后我做了收敛, 积分的高斯求积的阶数. 我选择了并使用域, 你可以改变以下。我有:
n direct rewritten
1 0.770878284949 0.770878284949
2 0.304480978430 0.304480978430
3 0.356922151260 0.356922151260
4 0.362576361509 0.362576361509
5 0.362316789057 0.362316789057
6 0.362314010897 0.362314010897
7 0.362314071949 0.362314071949
8 0.362314072182 0.362314072182
9 0.362314072179 0.362314072179
10 0.362314072179 0.362314072179
如您所见,对于两个积分都完全收敛到 12 位有效数字。
这是代码:
from scipy.integrate import fixed_quad
from scipy.special import jn
from numpy import exp, pi, sin, cos, array
def gauss(f, a, b, n):
"""Gauss quadrature"""
return fixed_quad(f, a, b, n=n)[0]
def f(x):
"""Function f(x) to integrate"""
return exp(-x) * x**2
xmax = 3.
print " n direct rewritten"
for n in range(1, 20):
def inner(theta_array):
return array([gauss(lambda x: f(x) * cos(x*sin(theta)), 0, xmax, n)
for theta in theta_array])
direct = gauss(lambda x: f(x) * jn(0, x), 0, xmax, n)
rewritten = gauss(inner, 0, pi, 20) / pi
print "%2d %.12f %.12f" % (n, direct, rewritten)
您可以自己玩这个,只需更改xmax
,可能您可能需要拆分间隔成元素,逐元素整合。您还可以更改功能f(x)
。确保始终收敛积分rewritten = gauss(inner, 0, pi, 20) / pi
,即从一些低阶开始并不断增加它,直到打印结果停止变化。