0

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

  • 1
    Is the indentation correctly reproduced? Then you are trying to solve the system before finishing to construct the `B` matrix. // You might also have to consider the idea that the condition number of the basis matrix becomes so bad with increasing size that the matrix is practically singular. – Lutz Lehmann Feb 12 '21 at 12:52
  • Thanks for the comment. I have edited my codes but another error arises. I will update my original post – Arnold Schwarzenegger Feb 12 '21 at 13:09
  • `a = fa + np.matmul(coef,B)`: You want the other monomial matrix there, with the powers of `x_list`. You might want to unify the variable names, `n_j, x_j, f_j` and `n_a,x_a,f_a` or similar. – Lutz Lehmann Feb 12 '21 at 13:29

0 Answers0