rbf.py

#!/usr/bin/python

import math

import numpy as np

from numpy import linalg as LA

# Function to approximate
def f(x):
    return math.exp(x * math.cos(3.0 * math.pi * x)) - 1.0

# Radial basis function (Gaussian with epsilon 3.05048)
def phi(r):
    return math.exp(-math.pow(3.05048 * r, 2.0))

# Discrete points
X = np.array([x / 14.0 for x in range(0, 15)])

# Data values (right hand side)
V = np.array([f(x) for x in X])

# Matrix containing interpolation values
# with phi(||xi - xj||) at each Mij
M = np.array([[phi(abs(x - y)) for x in X] for y in X])

print('X: %s' % X)
print('V: %s' % V)
print('M: %s' % M)

# Solve M*W=V
W = LA.solve(M, V)

print('W: %s' % W)

# Compute intepolated values
S = np.array([sum([w * phi(abs(x - y)) for w, y in zip(W, X)]) for x in X])

print('S: %s' % S)

# Maximum error
print('Error: %s' % max(np.array([abs(x - y) for x, y in zip(S, V)])))