4

我需要证明：

从 i=1 到 c sub i 绝对值平方的无穷大的和等于 1

烦人的是 c_i 等于函数 G 的积分。这是我的尝试。

import numpy as np
from scipy.integrate import quad

def G(x,n):
    P = (np.sqrt(735))*(np.sqrt(2))*np.sin(n*np.pi*x)*((x**3.0) - (11.0/7.0)*(x**2.0) + (4.0/7.0)*(x))
    return P

def Sum(x, n):
    i = 1
    S = 0
    I, err = quad(G, 0, 1, args=(i))
    while (i<n):
        S = S + I
        i = i + 1
    return S
x = np.linspace(0, 1, 250)    
print Sum(x, 200)

我遇到的问题是写总结部分。当我运行它时，我得到的数字越大，我给它的值就越多。如果选择 n 相当高（而不是无穷大），则可以显示总和如何趋于 1

4

2 回答 2

7

这个问题有很大的教学价值。

正如@alko 指出的那样，这个问题可以通过分析来解决。如果目的是证明总和等于一，那么应该分析地完成。

也许，这只是需要做的事情的一个简单版本，实际问题无法通过分析解决。在这种情况下，解决这样一个更简单的问题是很好的第一步。不幸的是，每当我们以数字方式解决问题时，我们都会引入新的问题集。

让我们按照问题中的建议和@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 的函数。它应该是一个。现在，当我们在 ipython 中运行它时，它会有点慢，并且第一次我们会收到一些长警告消息，但代码似乎会运行。此外，如果我们再次运行它（在同一个 ipython 会话中），我们将不会收到警告消息，所以我们可以忽略它，对吧？错了，但无论如何我们都会忽略它，因为这样做很常见。

如果我们看S它似乎在显示我们想要的东西，直到S[200]事情开始出错，价值开始增长！电脑怎么了？！什么都没有，我们忽略了另一个问题指标。 quad()返回误差估计值以及积分估计值。我们通常会忽略误差估计，但不应该这样做。如果我们绘制 S 和误差值，我们会发现以下内容。 G功能集成

所以我们看到是的，S 的值出现了可怕的错误，但quad()也告诉我们这会发生！事实上，我们忽略的警告也告诉了我们同样的事情。

我们如何理解并解决它？在这一点上，我将停止讲故事。如果盯着G(x,n)看不清楚，那么最好在n积分范围内绘制该函数。我们会发现它是一个剧烈振荡的函数，因此难以进行数值积分也就不足为奇了。一定会有更好的办法。

当然还有更好的方法。如果我们查看文档quad()并运行，quad_explain()我们可以了解权重函数。正弦是一种常见的权重函数，出现在积分中，因此有一些特殊的技术可以处理这种情况。因此，更好的方法如下（现在我们明白我为什么定义了g(x)：

def do_integral_weighted (N) :
    res = np.zeros(N)
    err = np.zeros(N)
    for n in xrange(1,N) :
        (res[n],err[n]) = integ.quad(g, 0, 1, weight='sin', wvar=n*np.pi)
    return (res,err)

(cw,errw) = do_integral_weighted(500)
Sw = np.cumsum(cw[1:]**2) # skip n=0

如图所示，我们发现它运行得更快并且更准确 g函数的集成 所以我们得到了一个稳定的答案，没有打印任何警告，并且有一个微小的集成错误。

我们从解决这个问题中学到了一些东西

  1. 可以分析解决的问题应该分析解决。
  2. 数学表达式的数值实现通常不像我们希望的那样简单。
  3. 不要忽略警告，也不要忽略我们正在使用的例程提供的错误信息。
  4. 数值技术不同于分析技术。numpy/scipy 提供了极其强大的工具。我们需要探索这些工具以充分利用它们的力量。
于 2013-10-25T14:55:36.237 回答
5

您的求和有误。在循环内部移动quad，并以这种方式计算平方和（而不是当前计算的 sum(c_i)）：

def Sum(x, n):
    S = 0
    for i in range(n):
        I, _err = quad(G, 0, 1, args=(i))
        S = S + I**2
    return S

顺便说一句，你怎么能证明这个总和的极限等于 1？它只能通过一些数学评估来显示，而不是计算。你可以说明这一点。

于 2013-10-24T22:30:36.530 回答