Code referenced in https://uli.rocks/p/polynomial-approximation/

sympy-exact-approx.py

# Run with ipython -i sympy-exact-approx.py so you can play with it afterwards

import sympy as sp

x = sp.Symbol('x')

# function to approximate, if this is exotic change inner() to compute
# integrals numerically instead of symbolically.
def f(x):
    # standard example
    return sp.sin(x)
    # exotic exmaple (can't be done symbolically). must change inner product
    # to compute integrals numerically.
    # return sp.exp(-(x*sp.cos(x))**2)

# inner product and bounds
a,b = -sp.pi, sp.pi
def inner(f, g):
    return sp.integrate(f*g, (x, a, b)) # symbolic
    # return sp.Integral(f*g, (x,a,b)).evalf(chop=True) # numeric

def normalize(v):
    return v / sp.sqrt(inner(v, v))

# degree of polynomial to use
deg = 5

B = [x**n for n in range(deg+1)] # basis
O = [normalize(B[0])]            # orthogonal basis

# orthogonalize B
O = [normalize(B[0])]
for i in range(1, len(B)):
    O.append(normalize(
        B[i] - sum(
          inner(B[i], O[j])*O[j] # <v_i, e_j>
          for j in range(0, i)
        )
    ))

# get coefficients
coeffs = [inner(f(x), O[j]) for j in range(len(O))]

# turn the coefficients into a sympy polynomial
poly = sum(O[j]*coeffs[j] for j in range(len(coeffs)))

# print the polynomial as latex!
print(sp.latex(poly))
# print the polynomial in numeric form
print(sp.latex(poly.evalf()))
# => 0.00564312 x^{5} - 0.1552714 x^{3} + 0.9878622 x


print(poly)
print('latex poly (for copying to desmos, etc)')
print(sp.latex(poly))

print('distance')
distance = sp.Integral((f(x)-poly)**2, (x,a,b)).evalf(chop=True)
print(distance)