在Python中编写无限和(具有整数的函数)

Question

我需要证明:

令人讨厌的是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

Answer 1

这个问题有很多教育价值.

正如@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的函数.它应该是一个.现在,当我们在ipython中运行它时,它会有点慢,我们第一次会得到一些长警告消息,但代码似乎会运行.此外,如果我们再次运行它(在同一个ipython会话中),我们将不会收到警告消息,所以我们可以忽略它,对吧？错了,但无论如何我们都会忽略它,因为这样做很常见.

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

所以我们看到是的,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

如图所示,我们发现运行速度更快,更准确 因此,我们得到一个稳定的答案,没有打印警告,以及一个微小的集成错误.

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

1. 分析可以解决的问题应该通过分析解决.
2. 数学表达式的数值实现通常不像我们希望的那样直截了当.
3. 不要忽略警告,也不要忽略我们正在使用的例程提供的错误信息.
4. 数值技术与分析技术不同.numpy/scipy提供了非常强大的工具.我们需要探索这些工具来充分利用它们的力量.

Answer 2

你的总结是错误的.移动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？它只能通过一些数学计算来显示,而不能通过计算来显示.你可以说明一下.