#===== Chebyshev interpolation for u(x) = sqrt(4+x) =====
import numpy as np
import scipy.integrate as sp
import matplotlib.pyplot as plt

n = 2
xk = np.array((-np.sqrt(3.)/2,0.,np.sqrt(3.)/2))  # n+1 Chebyshev nodes
uk = np.sqrt(4+xk)                                    # function values

#--- Newton divided difference
D = uk.copy()                       # array for the divided differences

for i in range(1,n+1):                  # start from 1 since D(0)=uu(0)
    for j in range(n,i-1,-1):                  # note the reverse order
        D[j]=(D[j] - D[j-1])/(xk[j] - xk[j-i])

print(D)

#--- compute the function, the Chebyshev interpolation and the error
x = np.linspace(-1,1,500)

u = np.sqrt(4+x)

I = D[0] + D[1]*(x-xk[0]) + D[2]*(x-xk[0])*(x-xk[1])

E = u - I

#--- plotting
plt.figure(1); plt.clf()
plt.plot(x,u,'k',label='$u(x)$')
plt.plot(x,I,'r-.',label='$I_2(x)$')
plt.xlabel('$x$')
plt.legend()

plt.figure(2); plt.clf()
plt.plot(x,E,'b')
plt.xlabel('$x$')
plt.ylabel('$E_2(x)$')
plt.axhline(linewidth=1,color='k')
plt.axvline(linewidth=1,color='k')

#--- Calculate c_3 in the Chebyshev least-squares approximation of u(x)
def integrand(y):
    return np.sqrt(4+np.cos(y))*np.cos(3*y)

ans, err = sp.quad(integrand,0,np.pi)

c3 = 2./np.pi*ans

print(c3)