Generating 1D recurrence coefficients for a given weight

orthopy contains four different algorithms that (in exact arithmetic) all do the same thing: They generate the recurrence coefficients for a 1D interval with a given weight function.

• Stieltjes.

  import orthopy
  import sympy
  
  N = 10
  
  alpha, beta = orthopy.tools.stieltjes(lambda t, ft: sympy.integrate(ft, (t, -1, 1)), N)
  print(alpha)
  print(beta)

  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
  [None, 1/3, 4/15, 9/35, 16/63, 25/99, 36/143, 49/195, 64/255, 81/323]
  

• Golub-Welsch. (Only works in floating-point arithmetic.) For this, you first need to compute the first 2*N moments of your measure

  integral(w(x) x^k dx)
  

  import quadpy
  
  N = 10
  
  moments = quadpy.tools.integrate(lambda x: [x ** (2 + k) for k in range(2 * N + 1)], -1, +1)
  moments = [float(m) for m in moments]  # convert to floats
  print(moments)
  print()
  
  alpha, beta = quadpy.tools.golub_welsch(moments)
  print(alpha)
  print(beta)

  [0.6666666666666666, 0.0, 0.4, 0.0, 0.2857142857142857, 0.0, 0.2222222222222222, 0.0, 0.18181818181818182, 
   0.0, 0.15384615384615385, 0.0, 0.13333333333333333, 0.0, 0.11764705882352941, 0.0, 0.10526315789473684,
   0.0, 0.09523809523809523, 0.0, 0.08695652173913043]
  
  [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
  [0.66666667 0.6        0.11428571 0.3968254  0.16161616 0.34265734
   0.18461538 0.31764706 0.19814241 0.30325815]
  

• Chebyshev. The moments can be used in exact arithmetic here.

  import quadpy
  
  N = 10
  
  moments = quadpy.tools.integrate(lambda x: [x ** (2 + k) for k in range(2 * N)], -1, +1)
  print(moments)
  print()
  
  alpha, beta = quadpy.tools.chebyshev(moments)
  print(alpha)
  print(beta)

  [2/3 0 2/5 0 2/7 0 2/9 0 2/11 0 2/13 0 2/15 0 2/17 0 2/19 0 2/21 0]
  
  [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
  [2/3, 3/5, 4/35, 25/63, 16/99, 49/143, 12/65, 27/85, 64/323, 121/399]
  

• Modified Chebyshev. For this, you first need to compute the first 2*N modified moments of your measure

  integral(w(x) p_k(x) dx)
  

  with a particular set of orthogonal polynomials p_k. A common choice are the Legendre polynomials, here taken from orthopy.

  import orthopy
  import quadpy
  
  N = 10
  
  
  def leg_polys(x):
      return orthopy.line_segment.tree_legendre(x, 2 * N - 1, "monic", symbolic=True)
  
  
  moments = quadpy.tools.integrate(
      lambda x: [x ** 2 * leg_poly for leg_poly in leg_polys(x)], -1, +1
  )
  print(moments)
  print()
  
  _, _, a, b = orthopy.line_segment.recurrence_coefficients.legendre(
      2 * N, "monic", symbolic=True
  )
  alpha, beta = quadpy.tools.chebyshev_modified(moments, a, b)
  print(alpha)
  print(beta)

  [2/3 0 8/45 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
  
  [0 0 0 0 0 0 0 0 0 0]
  [2/3 3/5 4/35 25/63 16/99 49/143 12/65 27/85 64/323 121/399]
  

Note: If you perform the computation in floating-point arithmetic for example because you need to compute the scheme fast), it makes sense to think about the numerical implications here. That's because the map from the moments to the recurrence coefficients can be very ill-conditioned, meaning that small round-off errors can lead to very wrong coefficients. For further computation, it's numerically beneficial if the moments are either 0 or in the same order of magnitude. That is why the flexibility provided by the modified Chebyshev scheme can be important. (Note how almost all moments are 0 there.)