-- An arbitrary function can be approximated by a weighted sum of Chebyshev polynomials.
-- A Chebyshev polynomial of degree n is defined as:
-- T(x)[n] = 1 when n=0,
--           x when n=1,
--           2*x*T(x)[n-1] - T(x)[n-2]
-- or the explicit formula:

T(x)[n] = cos(n*acos(x))

Xmin=-1; Xmax=1; Title="Chebyshev polynomials"
plot T(X)[0:6]  -- show the first few polynomials

-- The weighting coefficients can be calculated using the values of the arbitrary function at the zeros of T(x). T(x)[n] has n zeros between -1 and +1 located at:

xz[k] = cos(pi*Radians*(k-.5)/n) dim[n]

-- the N coefficients:
c[j] = (2/n)*sum(f(xz[k])*T(xz[k])[j-1],k,1,n) dim[n]
cc:=c:;  -- save the coefficients

-- the approximation formula:
approx(cc,x) = sum(cc[k]*T(x)[k-1],k,1,count(cc)) - cc[1]/2

-- Try an example function. Smooth functions are easy. Try something a bit tougher.

f(x) = x when x<-.5, -.5 when x<.5, -x

n=15;  -- number of coefficients to use

newaxis
plot f(X)
plot {xz,f(xz)}
plot approx(cc,X)

-- Once the coefficients are known they can be used to derive coefficients to approximate the integral.

ci[i] = (cc[i-1]-cc[i+1])/(2*(i-1)) when i>1,
          ci0 dim[n-1]
ci0 = 0   -- arbitrary constant of integration

newaxis 
plot approx(ci,X)

-- They can also be used to find coefficients to approximate the the derivative.

cd:=0[1:n]:
cd[n-1]:=2*(n-1)*cc[n]:
j:=n-2:
(cd[j]:=cd[j+2] + 2*j*cc[j+1], j:=j-1) while j>0:

plot approx(cd,X)

