In [8]:
from __future__ import division
import numpy as np
import scipy as sp
import scipy.special as ss
import matplotlib.pyplot as pt
import numpy.linalg as la

%pylab inline

Populating the interactive namespace from numpy and matplotlib

WARNING: pylab import has clobbered these variables: ['f']
`%pylab --no-import-all` prevents importing * from pylab and numpy
In [21]:
nelements = 5
nnodes = 3
mesh = np.linspace(-1, 1, nelements+1, endpoint=True)
gauss_nodes = ss.legendre(nnodes).weights[:, 0]*0.5 + 0.5

widths = np.diff(mesh)
nodes = mesh[:-1, np.newaxis] + widths[:, np.newaxis] * gauss_nodes
In [22]:
def f(x):
    return np.abs(x-0.123812378)

pt.plot(nodes.flat, f(nodes).flat)
Out[22]:
[<matplotlib.lines.Line2D at 0x107b3dc10>]
In [23]:
nmany_nodes = 32

many_gauss_nodes = ss.legendre(nmany_nodes).weights[:, 0]*0.5 + 0.5
many_nodes = mesh[:-1, np.newaxis] + widths[:, np.newaxis] * many_gauss_nodes

def legendre_vdm(nodes, nmodes):
    result = np.empty((len(nodes), nmodes))
    for i in xrange(nmodes):
        result[:, i] = ss.eval_legendre(i, nodes)
    return result

vdm = legendre_vdm(gauss_nodes, nnodes)
many_vdm = legendre_vdm(many_gauss_nodes, nnodes)
zero_pad = np.zeros((nmany_nodes, nnodes))
zero_pad[:nnodes, :nnodes] = np.eye(nnodes)
upterpolate = np.dot(many_vdm, la.inv(vdm))
In [24]:
fnodes = f(nodes)
fmany_nodes = np.dot(upterpolate, fnodes.T).T
pt.plot(many_nodes.flat, fmany_nodes.flat)
Out[24]:
[<matplotlib.lines.Line2D at 0x107bcc410>]
In [ ]: