Writing an infinite sum (of a function that has an integral) in Python

Question

I am required to show that:

The annoying thing is that c_i is equal to the integral of the function G. Here is my attempt.

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)

The problem I am having is writing the summation part. When I run this I get a bigger number the more values I give it. If n is chosen to be quite high (instead of infinity) it can be shown how the sum tends to 1

Answer 1

This problem has much pedagogic value.

As @alko points out, this problem can be solved analytically. If the intention is to prove that the sum is equal to one then it should be done analytically.

Perhaps, this is just a simpler version of something that needs to be done and the actual problem is not solvable analytically. In that case, solving a simpler problem such as this one is a good first step. Unfortunately, whenever we solve something numerically we introduce new sets of problems.

Let us proceed as suggested in the question and the correction given by @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

(The reason why I define two "G" functions will be apparent below.)

Once run the array S will contain the desired sum as a function of n. It should be one. Now, when we run this within ipython it will be somewhat slow and the first time we will get some long warning message, but the code will seem to run. Further, if we run it again (in the same ipython session) we will not get a warning message so we can just ignore that, right? Wrong but we will ignore it anyway since it is common to do so.

If we look at S it seems to be showing what we want until about S[200] where things start to go wrong, the value starts to grow! What is wrong with the computer?! Nothing, and we ignored another indicator of a problem. quad() returns an error estimate along with an estimate of the integral. We typically ignore the error estimate but shouldn't. If we plot the S and error values we find the following.

So we see that yes, the value of S goes horribly wrong, but quad() was also telling us this would happen! In fact, the warning we ignored also told us the same thing.

How do we understand this and fix it? At this point I will stop with the story telling. If staring at the G(x,n) doesn't make it clear then it would be good to plot that function for a large n over the range of integration. We will find it is a wildly oscillating function so it isn't surprising that it is hard to integrate numerically. There must be a better way.

Of course there is a better way. If we look at the documentation for quad() and run quad_explain() we can learn about weight functions. The sine is a common weight function that appears in integrals so there are special techniques for handling this case. Thus, a better approach is as follows (now we see why I defined 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

This we find runs much faster and is much more accurate as shown in the plot So we get a stable answer, no warnings printed, and a tiny integration error.

We learn a few things from working this problem

1. Problems that can be solved analytically should be solved analytically.
2. The numerical implementation of a mathematical expression is often not as straight forward as we may have hoped.
3. Do not ignore warnings and do not ignore the error information provided by the routines we are using.
4. Numerical techniques are different than analytic techniques. numpy/scipy provide extremely powerful tools. We need to explore these tools to fully exploit their power.

Answer 2

You have errror in your summing. Move quad inside loop, and evaluate sum of squares (instead of currently evaluated sum(c_i)) this way:

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

By the way, how can you SHOW that limit of this sum equals 1? It can be only shown by some math evaluations, not computations. You can just illustrate that.