这个问题很有教育价值。在
正如@alko指出的,这个问题可以通过分析来解决。如果目的是证明总和等于1,则应进行分析。在
也许,这只是需要做的事情的一个更简单的版本,而实际的问题是不可解析地解决的。在这种情况下,解决一个像这样的简单问题是一个好的第一步。不幸的是,每当我们用数值方法解决某个问题时,我们就会引入一系列新的问题。在
让我们按照问题中的建议和@alko给出的更正进行。在import numpy as np
import scipy.integrate as integ
def g(x) :
return np.sqrt(1470) * (x**3-11./7*x**2+4./7*x)
def G(x,n) :
return np.sin(n*np.pi*x) * g(x)
def do_integral (N) :
res = np.zeros(N)
err = np.zeros(N)
for n in xrange(1,N) :
(res[n],err[n]) = integ.quad(G, 0, 1, args=(n,))
return (res,err)
(c,err) = do_integral(500)
S = np.cumsum(c[1:]**2) # skip n=0
(我定义两个“G”函数的原因如下所示。)
一旦运行数组S将包含所需的和作为n的函数。它应该是1。现在,当我们在ipython中运行它时,它会有点慢,第一次我们会收到一些很长的警告消息,但是代码似乎会运行。此外,如果我们再次运行它(在同一个ipython会话中),我们将不会收到警告消息,所以我们可以忽略它,对吗?错误但我们无论如何都会忽略它,因为这样做很常见。在
如果我们看一下S,它似乎显示了我们想要什么,直到S[200]事情开始出错时,值开始增长!电脑怎么了?!没什么,我们忽略了另一个问题的迹象。quad()返回一个错误估计值和一个积分估计值。我们通常忽略误差估计,但不应该忽略。
所以我们看到,是的,S的值出了严重的错误,但是quad()也告诉我们这会发生!事实上,我们忽视的警告也告诉了我们同样的事情。在
我们如何理解并修复它?在这一点上,我将停止讲故事。如果盯着G(x,n)看不清楚,那么最好在积分范围内绘制一个大的n函数。我们会发现它是一个剧烈振荡的函数,所以很难用数值积分也就不足为奇了。一定有更好的办法。在
当然还有更好的办法。如果我们查看quad()的文档并运行quad_explain(),我们可以了解权重函数。正弦函数是一种常见的权函数,出现在积分中,所以有特殊的技术来处理这种情况。因此,一个更好的方法如下(现在我们来看看为什么我定义g(x)):
^{pr2}$
我们发现这一点运行得更快,也更准确,如图所示
所以我们得到了一个稳定的答案,没有打印警告,还有一个微小的集成错误。在
我们从解决这个问题中学到了一些东西可以分析解决的问题应该分析解决。在
数学表达式的数值实现往往不像我们希望的那样直接。在
不要忽略警告,也不要忽略我们正在使用的例程提供的错误信息。在
数值技术不同于分析技术。numpy/scipy提供了非常强大的工具。我们需要探索这些工具来充分利用它们的力量。在