In [127]:
import matplotlib.pyplot as plt
import numpy as np
import seaborn 
seaborn.set(style='whitegrid',rc={'figure.figsize':(9,5)})

(a) Fit the data with a polynomial with 5, 10, and 15 terms, using a pseudo-inverse of the Vandermonde matrix.

In [128]:
x_data = np.arange(-10,11)
y_data = np.zeros(len(x_data))

for i in range(len(x_data)):
    if x_data[i]>0: y_data[i]=1
In [129]:
VSolvePlot(x_data,y_data,5)
VSolvePlot(x_data,y_data,10)
VSolvePlot(x_data,y_data,15)

5 terms

10 terms

15 terms
In [130]:
def VSolvePlot(x_,y_,n_):
    n_plot = 10000
    x_plot = np.linspace(-10,10,n_plot)
    y_plot = np.zeros(n_plot)
    A = np.empty([len(x_),n_])
    for i in range(len(x_)):
        for j in range(n_):
            A[i][j] = x_[i]**j
    A_pseudo = np.linalg.pinv(A, rcond=1e-15,)
    v_terms = np.dot(A_pseudo,y_)
    for i in range(n_plot):
        y_plot[i] = sum(v_terms[j]*(x_plot[i]**j) for j in range(n_))
    
    print(n_,"terms")
    fig,ax = plt.subplots()
    plt.ylim(-0.2,1.2)
    plt.plot(x_plot,y_plot,'black',lw=1)
    ax = plt.scatter(x_,y_,40,'blue')
    plt.show()

(b)

In [135]:
RBFSolve(x_data,y_data,5)
RBFSolve(x_data,y_data,10)
RBFSolve(x_data,y_data,15)

5 terms

10 terms

15 terms
Out[135]:
(array([ 0.00273591, -0.00808859,  0.01471804, -0.02715188,  0.05045257,
        -0.09385279,  0.1725193 , -0.21704709,  0.11759003,  0.0246744 ,
        -0.05932695,  0.03316382, -0.01318932,  0.00143943,  0.00106801]),
 array([-10.        ,  -8.66666667,  -7.33333333,  -6.        ,
         -4.66666667,  -3.33333333,  -2.        ,  -0.66666667,
          0.66666667,   2.        ,   3.33333333,   4.66666667,
          6.        ,   7.33333333,   8.66666667]))
In [137]:
def RBFFunc(x, a, c):
    n_ = len(a)
    y = sum(a[j]*abs(x-c[j])**3 for j in range(n_))
    return y
In [138]:
def RBFSolve(x_,y_,n_):
    c = np.linspace(-10,10,n_,endpoint=False)
    n_plot = 10000
    x_plot = np.linspace(-10,10,n_plot)
    y_plot = np.zeros(n_plot)
    A = np.empty([len(x_),n_])
    for i in range(len(x_)):
        for j in range(n_):
            A[i][j] = abs(x_[i]-c[j])**3
    A_pseudo = np.linalg.pinv(A, rcond=1e-17,)
    v_coef = np.dot(A_pseudo,y_)
    for i in range(n_plot):
        y_plot[i] = RBFFunc(x_plot[i], v_coef, c)
    print(n_,"terms")
    
    fig,ax = plt.subplots()
    plt.ylim(-0.2,1.2)
    plt.plot(x_plot,y_plot,'black',lw=1)
    ax = plt.scatter(x_,y_,40,'blue')
    plt.show()
    
    return v_coef, c

(c)

In [153]:
x_outofsamp = np.arange(-10.5,10.5)
RBFError(x_outofsamp,5)
RBFError(x_outofsamp,10)
RBFError(x_outofsamp,15)

5 terms

error = 0.1591072612797621
10 terms

error = 0.1319499010665493
15 terms

error = 0.10155710585145762
In [151]:
def RBFError(x_,n_):
    v_coef, c = RBFSolve(x_data,y_data,n_)
    y_fit = np.empty(len(x_))
    for i in range(len(x_)):
        y_fit[i] = RBFFunc(x_[i], v_coef, c)
    y_exact = BaseFunc(x_)
    error = np.sqrt(sum((y_fit[i]-y_exact[i])**2 for i in range(len(x_)))/len(x_))
    print("error =", error)
In [149]:
def BaseFunc(x_):
    y_ = np.zeros(len(x_))
    for i in range(len(x_)):
        if x_[i]>0: y_[i]=1
    return y_