title: python: integrating a piecewise function
body: <p>I want to integrate a piecewise a defined function that is multiplied by the Legendre polynomials.
  Unfortunately, I can't find how to use the nth Legendre polynomial of x in the <a href="http://docs.scipy.org/doc/numpy/reference/routines.polynomials.legendre.html" rel="nofollow">documentation</a>.  I want to integrate each Legendre polynomial of x when <code>n = 1,..., 50</code> so I have set <code>n = np.arange(1, 51, 1)</code>.</p>
  
  <pre><code>import numpy as np
  import pylab
  from scipy import integrate
  
  n = np.arange(1, 51, 1)                                                   
  
  
  def f(x):
      if 0 &lt;= x &lt;= 1:
          return 1
      if -1 &lt;= x &lt;= 0:
          return -1
  </code></pre>
  
  <p>I suppose I need to define another function let's say <code>u(x)</code>.</p>
  
  <pre><code>c = []
  
  
  def u(x):
      c.append((2. * n + 1) / 2. * integrate.quad(f(x) * insert Legendre polynomials here, -1., 1.)) 
      return sum(c * Legendre poly, for nn in range(1, 51)) 
  </code></pre>
  
  <p>So I would then return some <code>u(x)</code> with the first 50 terms expanding my piecewise function by Legendre polynomials.</p>
  
  <p><strong>Edit 1:</strong></p>
  
  <p>If this can't be done, I could use Rodrigues's Formula to compute the nth Legendre polynomial.  However, I couldn't find anything useful when I was looking for computing nth derivatives in Python.</p>
  
  <pre><code>P_n(x) = \frac{1}{2^n n!}\frac{d^n}{dx^n}(x^2 - 1)^n
  </code></pre>
  
  <p>So this is an option if someone knows how to implement such a scheme in Python.</p>
  
  <p><strong>Edit 2:</strong></p>
  
  <p>Using Saullo Castro's answer, I have:</p>
  
  <pre><code>import numpy as np
  from scipy.integrate import quad
  
  def f(x, coef):
      global p
      p = np.polynomial.legendre.Legendre(coef=coef)
      if 0 &lt;= x &lt;= 1:
          return 1*p(x)
      if -1 &lt;= x &lt;= 0:
          return -1*p(x)
  
  c = []
  for n in range(1, 51):
      c.append((2. * n + 1.) / 2. * quad(f, -1, 1, args=range(1,n+1))[0])
  
  def g(x)
      return sum(c * p(x) for n in range(1, 51))
  </code></pre>
  
  <p>However, if I print <code>c</code>, the values are wrong.  The values should be <code>1.5, 0, -7/8, 0, ...</code></p>
  
  <p>Also, when I plot <code>g</code>, I would like to do <code>x = np.linspace(-1, 1, 500000)</code> so the plot is detailed but <code>c</code> is only 50.  How can this be achieved?</p>