Skip to content

Instantly share code, notes, and snippets.

Created June 16, 2023 19:42
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
What would you like to do?
# Run with ipython -i 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)):
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 the polynomial in numeric form
# => 0.00564312 x^{5} - 0.1552714 x^{3} + 0.9878622 x
print('latex poly (for copying to desmos, etc)')
distance = sp.Integral((f(x)-poly)**2, (x,a,b)).evalf(chop=True)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment