text: #!/usr/bin/env python
  
  """
  This example illustrates the Gibbs phenomenon.
  
  It also calculates the Wilbraham-Gibbs constant by two approaches:
  
  1) calculating the fourier series of the step function and determining the
  first maximum.
  2) evaluating the integral for si(pi).
  
  See:
   * https://en.wikipedia.org/wiki/Gibbs_phenomena
  """
  
  from sympy import var, sqrt, integrate, conjugate, seterr, Abs, pprint, I, pi,\
      sin, cos, sign, lambdify, Integral, S
  
  x = var("x", real=True)
  
  
  def l2_norm(f, lim):
      """
      Calculates L2 norm of the function "f", over the domain lim=(x, a, b).
  
      x ...... the independent variable in f over which to integrate
      a, b ... the limits of the interval
  
      Examples
      ========
  
      >>> from sympy import Symbol
      >>> from gibbs_phenomenon import l2_norm
      >>> x = Symbol('x', real=True)
      >>> l2_norm(1, (x, -1, 1))
      sqrt(2)
      >>> l2_norm(x, (x, -1, 1))
      sqrt(6)/3
  
      """
      return sqrt(integrate(Abs(f)**2, lim))
  
  
  def l2_inner_product(a, b, lim):
      """
      Calculates the L2 inner product (a, b) over the domain lim.
      """
      return integrate(conjugate(a)*b, lim)
  
  
  def l2_projection(f, basis, lim):
      """
      L2 projects the function f on the basis over the domain lim.
      """
      r = 0
      for b in basis:
          r += l2_inner_product(f, b, lim) * b
      return r
  
  
  def l2_gram_schmidt(list, lim):
      """
      Orthonormalizes the "list" of functions using the Gram-Schmidt process.
  
      Examples
      ========
  
      >>> from sympy import Symbol
      >>> from gibbs_phenomenon import l2_gram_schmidt
  
      >>> x = Symbol('x', real=True)    # perform computations over reals to save time
      >>> l2_gram_schmidt([1, x, x**2], (x, -1, 1))
      [sqrt(2)/2, sqrt(6)*x/2, 3*sqrt(10)*(x**2 - 1/3)/4]
  
      """
      r = []
      for a in list:
          if r == []:
              v = a
          else:
              v = a - l2_projection(a, r, lim)
          v_norm = l2_norm(v, lim)
          if v_norm == 0:
              raise ValueError("The sequence is not linearly independent.")
          r.append(v/v_norm)
      return r
  
  
  def integ(f):
      return integrate(f, (x, -pi, 0)) + integrate(-f, (x, 0, pi))
  
  
  def series(L):
      """
      Normalizes the series.
      """
      r = 0
      for b in L:
          r += integ(b)*b
      return r
  
  
  def msolve(f, x):
      """
      Finds the first root of f(x) to the left of 0.
  
      The x0 and dx below are tailored to get the correct result for our
      particular function --- the general solver often overshoots the first
      solution.
      """
      f = lambdify(x, f)
      x0 = -0.001
      dx = 0.001
      while f(x0 - dx) * f(x0) > 0:
          x0 = x0 - dx
      x_max = x0 - dx
      x_min = x0
      assert f(x_max) > 0
      assert f(x_min) < 0
      for n in range(100):
          x0 = (x_max + x_min)/2
          if f(x0) > 0:
              x_max = x0
          else:
              x_min = x0
      return x0
  
  
  def main():
      L = [1]
      for i in range(1, 100):
          L.append(cos(i*x))
          L.append(sin(i*x))
      # next 2 lines equivalent to L = l2_gram_schmidt(L, (x, -pi, pi)), but faster:
      L[0] /= sqrt(2)
      L = [f/sqrt(pi) for f in L]
  
      f = series(L)
      print("Fourier series of the step function")
      pprint(f)
      x0 = msolve(f.diff(x), x)
  
      print("x-value of the maximum:", x0)
      max = f.subs(x, x0).evalf()
      print("y-value of the maximum:", max)
      g = max*pi/2
      print("Wilbraham-Gibbs constant        :", g.evalf())
      print("Wilbraham-Gibbs constant (exact):", \
          Integral(sin(x)/x, (x, 0, pi)).evalf())
  
  if __name__ == "__main__":
      main()
  