import numpy as np
import numpy.linalg as la
import matplotlib.pyplot as plt
Here's the half-circle function:
def f(x):
return np.sqrt(1-x**2)
We'll only look at it on $[0,1]$:
x_plot = np.linspace(0, 1, 1000)
plt.plot(x_plot, f(x_plot))
degree = 4
nodes = np.linspace(0, 1, degree+1)
nodes
Now build the Vandermonde matrix:
V = np.array([
nodes**i
for i in range(degree+1)
]).T
And find the coefficients as coeffs
:
coeffs = la.solve(V, f(nodes))
Now evaluate the interpolant at x_plot
interp_plot = 0*x_plot
for i in range(degree+1):
interp_plot += coeffs[i] * x_plot**i
plt.plot(x_plot, interp_plot)
plt.plot(x_plot, f(x_plot))
plt.plot(nodes, f(nodes), "o")
Now integrate the interpolant:
# clear
integral = 0
for i in range(degree+1):
integral += coeffs[i] * 1/(i+1) * (1**(i+1) - 0**(i+1))
my_pi = 4*integral
my_pi