Here is my attempt to perform polynomial interpolation with monomial basis B = {1,x,x^2,x^3,...x^n}. I keep getting errors while solving for the coefficients from the linear system cB = f
Here are my codes:
import numpy as np
import matplotlib.pyplot as plt
nvec = np.array([3, 10, 20, 50]) # n intervals as given in the problem
xL = 1; xR = 4; #left and right bounds
def f(x):
return np.sin(np.cos(3*x))
error = np.zeros(len(nvec)) # a row of zeros with length=4
plt.figure()
for j in range(len(nvec)):
n = nvec[j]
npoints = np.arange(0,n+0.1,1) #Chebyshev points
xj=1/2*(xL+xR) + 1/2*(xL-xR)*np.cos(npoints*np.pi/n);
#Chebyshev point distribution
fj = f(xj) #value of true function at each xj
n_a = 1001 #Number of points to plot approximation
x_list = np.linspace(xL, xR, n_a)
fa = np.zeros(n_a)
B = np.zeros([n+1,n+1]) #empty matrix for monomial basis
for i in range(0,n+1):
B[::,i] = np.power(xj.T,i) #raise the elements in each column evaluated
# at xj to the power of i -- {1,x,x^2,x^3,etc}
#print(B.shape)
coef = np.linalg.solve(B,fj)
coef = np.flip(coef, axis=0)
print(coef.shape)
print(B.shape)
fa = fa + np.matmul(coef,B)
error[j] = np.linalg.norm(f(x_list) - fa[x_list]);
ipl = '22'+str(j+1)
plt.subplot(ipl)
plt.ylim(-1,1)
plt.xlim(1,4)
plt.plot(x_list, f(x_list), linestyle='-', color='g', label='$f (x)$')
plt.plot(x_list, fa, linestyle='--', color='r', label='$fa (x)$')
plt.tick_params(axis="both",direction="in")
plt.tick_params(right=True)
plt.scatter(xj, f(xj), marker='o', s=10, color='k', label='$f (x_j)$' )
plt.title('$n=%i$'%(nvec[j]))
#plt.axis(’equal’)
if j<2:
plt.xticks([])
if j==2 or j==3:
plt.xlabel('$x$')
if j==0:
plt.legend()
plt.subplots_adjust( wspace=0.3, hspace=0.25)
plt.figure(figsize=(10,12))
plt.rcParams['font.size']=13
plt.gca().set_yscale('log')
plt.title('Interpolation error')
plt.xlabel('$n$'); plt.ylabel('$||f(x)-f_a(x)||$')
plt.xlim(0, nvec[len(nvec)-1])
plt.tick_params(right=True)
plt.scatter(nvec, error, s=100, color='y')
plt.tick_params(axis="both",direction="in")
plt.tick_params(right=True, top=True)
Here is the error message :
a = fa + np.matmul(coef,B)
ValueError: operands could not be broadcast together with shapes (1001,) (4,)
<Figure size 432x288 with 0 Axes>
I can't figure out how to fix the error here/ I have set up a N+1 by N+1 matrix B where each column corresponds to an elements in the monomial basis {1,x,x^2,.....} . I have the exact values of a function f at each xj written as fj, and I just need to solve for the coefficients. It is weird to me that monomial matrix B cannot be multiplied by the coefficient matrix coef