Galerkin method¶

The cubic B-spline basis¶

In [7]:
import sympy as sp
sp.init_printing()
In [8]:
import numpy as np
import matplotlib.pyplot as plt
In [9]:
a = np.array([1,2,3,4,5])
b = np.array([True,True,False,False,True])
a[b] = 999
a
Out[9]:
array([999, 999,   3,   4, 999])
In [10]:
def cubic_bspline(x):
    f = np.zeros_like(x)
    subinterval = (x>=-2) & (x<-1); subx = x[subinterval]+2; f[subinterval] = subx**3
    subinterval = (x>=-1) & (x< 0); subx = x[subinterval]+1; f[subinterval] = 1 + 3*subx + 3*subx**2 - 3*subx**3
    subinterval = (x>= 0) & (x< 1); subx = 1-x[subinterval]; f[subinterval] = 1 + 3*subx + 3*subx**2 - 3*subx**3
    subinterval = (x>= 1) & (x< 2); subx = 2-x[subinterval]; f[subinterval] = subx**3
    return f

x = np.linspace(-3,3,200)
plt.plot(x,cubic_bspline(x));
In [14]:
# Example of a function that creates and returns a function
def make_adder(n):
    def adder(x):
        return x+n
    return adder
f = make_adder(3)
g = make_adder(100)
f(5),g(5)
Out[14]:
$\displaystyle \left( 8, \ 105\right)$
In [17]:
def make_phi(template,h,xi):
    def phi(x):
        return template( (x-xi)/h )
    return phi

def make_phi(template,h,xis,cis):  # allow linear combinations to handle ends of domain
    def phi(x):
        return np.sum([ci*template( (x-xi)/h ) for ci,xi in zip(cis,xis)],axis=0) 
    return phi

gridx = np.linspace(0,1,10)
h = gridx[1] - gridx[0]
phis = [make_phi(cubic_bspline,h,gridx[:2],[-1/4,1])] + \
       [make_phi(cubic_bspline,h,[xi],[1]) for xi in gridx[2:-2]] + \
       [make_phi(cubic_bspline,h,gridx[-2:],[1,-1/4])] 
x = np.linspace(0,1,500)
[plt.plot(x, phi(x)) for phi in phis];

Examples of linear combinations of cubic B-splines¶

In [18]:
for example in range(1):
    plt.plot(x, np.sum([np.random.rand()*phi(x) for phi in phis],axis=0) )
In [19]:
plt.plot(x, np.sum([i*phi(x) for i,phi in enumerate(phis[3:12])],axis=0) );

Approximation of some given smooth functions by cubic B-splines¶

In [22]:
# do a cubic b-spline approximation to a given function
# by making the difference between the approximation and the function orthogonal to the basis functions

from scipy import integrate
np.set_printoptions(linewidth=300)

def f1(x): return np.sin(2*np.pi*x) #x*(1-x)

n = len(phis)
A = np.zeros((n,n))
b = np.empty(n)
for i in range(n):
    for j in range(n):
        F = lambda x:phis[i](x)*phis[j](x)
        A[i,j] = integrate.quad(F, 0, 1)[0]  # make sure adaptive quadrature is picking up localized blips
        if i!=j: A[j,i] = A[i,j]
    F = lambda x: f1(x)*phis[i](x)
    b[i] = integrate.quad(F,0,1)[0]
A            
b
c = np.linalg.solve(A,b)
plt.figure(figsize=(5,8))
plt.subplot(2,1,1)
plt.plot(x,f1(x),label='function to be approximated')
approx = np.sum([ci*phi(x) for ci,phi in zip(c,phis)],axis=0)
plt.plot(x, approx,label='cubic spline approximation')
plt.legend()
plt.subplot(2,1,2)
plt.plot(x,approx-f1(x),label='error')
plt.legend();
In [24]:
# do a cubic b-spline approximation to a given function
# by making the difference between the approximation and the function orthogonal to the basis functions

from scipy import integrate
np.set_printoptions(linewidth=300)

def f1(x): return  x**2*(1-x)

n = len(phis)
A = np.zeros((n,n))
b = np.empty(n)
for i in range(n):
    for j in range(n):
        F = lambda x:phis[i](x)*phis[j](x)
        A[i,j] = integrate.quad(F, 0, 1)[0]  # make sure adaptive quadrature is picking up localized blips
        if i!=j: A[j,i] = A[i,j]
    F = lambda x: f1(x)*phis[i](x)
    b[i] = integrate.quad(F,0,1)[0]
A            
b
c = np.linalg.solve(A,b)
plt.figure(figsize=(5,8))
plt.subplot(2,1,1)
plt.plot(x,f1(x),label='function to be approximated')
approx = np.sum([ci*phi(x) for ci,phi in zip(c,phis)],axis=0)
plt.plot(x, approx,label='cubic spline approximation')
plt.legend()
plt.subplot(2,1,2)
plt.plot(x,approx-f1(x),label='error')
plt.legend();
In [ ]: