In this lecture we'll try to construct an L^2 projection using some of the objects defined in the project ePICURE: http://github.com/luca-heltai/ePICURE
In particular we will be using some specializations of the VectorSpace class. The UniformLagrangeVectorSpace with N degrees of freedom.
Let's start by loading all required packages and libraries.
In [ ]:
%pylab inline
import interfaces
reload(interfaces)
from interfaces import *
Now let's define a small helper function which will help us visualize what's going on.
In [ ]:
def plot_all(vs, i, span=1025):
"""Plot all basis of the VectorSpace vs in figure i using span points (by default: 1025)"""
x = linspace(vs.cells[0], vs.cells[-1], span)
figure(i)
for i in range(vs.n_dofs):
plot (x, vs.basis(i)(x))
title(type(vs).__name__)
show()
In [4]:
N = 3
vs = UniformLagrangeVectorSpace(N)
plot_all(vs, 1)
What else do we know about this VectorSpace? Let's try its print_info() method:
In [ ]:
vs.print_info()
As we can see, there are three degrees of freedom, 1 cell (with cell boundaries [0,1]) and in the only cell we have, we get all indices of the nonzero basis functions.
What if we wanted to transform this basis function to a different interval? Say [.3, .55]?
There is a special VectorSpace, called AffineVectorSpace, which transforms arbitrary VectorSpaces affinely, to the interval given as argument: let's try this...
In [6]:
va = AffineVectorSpace(vs, .3, .55)
plot_all(va, 2)
Let's see what has changed:
In [ ]:
va.print_info()
In [8]:
newq = linspace(0,1,3)
vi = IteratedVectorSpace(vs, newq)
plot_all(vi, 2)
In [ ]:
vi.print_info()
As we can see, now things have changed. We used newq as a definition of subintervals, and we asked IteratedVectorSpace to repeat the VectorSpace passed as argument as many times as there are intervals in newq. What happens at the elements boundaries? The UniformLagrangeVectorSpace has an additional argument at construction time, which determines how many dofs are shared across cell boundaries. Let's try to set that argument to False, and repeat the construction:
In [ ]:
vs_dg = UniformLagrangeVectorSpace(3, False)
vs_dg.print_info()
Now you can see that the number of shared dofs across boundaries have changed. Let's analyize how this will affect the IteratedVectorSpace:
In [11]:
vi_dg = IteratedVectorSpace(vs_dg, newq)
plot_all(vi_dg, 3)
They seem almost identical. Look closely, and you'll notice that now there are more colors, and that some functions which were continuous before, now are no longer continuous...
Let's see how many degrees of freedom we have:
In [ ]:
vi_dg.print_info()
You can see now that there are no overlaps between cells. All basis functions are local. Let's try to construct an interpolation matrix and an L^2 projection matrix, and see how these changes will affect the solution...
We have already seen how to construct a least square matrix. I have implemented this in the utilities folder, in the least_square.py module.
Lets load it. Notice that I put a reload command right after import utilities.least_square. This allows me to make changes in the given file, and let this notebook notice the changes...
In [ ]:
import utilities.matrices
reload(utilities.matrices)
from utilities.matrices import *
# This will create the matrix M_ij = v_j(x_i), on 50 interpolation points
x = linspace(0,1,50)
M = InterpolationMatrix(vi, x)
type(M)
Notice that we can plot the matrix, if we want... it should look similar to the plots above, just with less points... Notice that the matrix above is a Column oriented matrix, because it was faster to produce than a row oriented matrix...
In [ ]:
# This will create the matrix M_ij = v_j(x_i), on 50 interpolation points
x = linspace(0,1,50)
M = InterpolationMatrix(vi, x)
type(M)
In [15]:
plot(x, M)
Out[15]:
There are a lot of zero entries in the M matrix. We can make it sparse:
In [ ]:
from scipy.sparse import csc_matrix
Ms = csc_matrix(M)
Now let's try to solve a least square problem with a given curve... The curve is a circle.
In [17]:
from numpy.linalg import lstsq
fx = lambda x: sin(2*pi*x)
fy = lambda x: cos(2*pi*x)
g = asmatrix(array([fx(x), fy(x)]))
c = lstsq(M, g.T)
# The first argument of c contains a matrix. Transform it to an array,
# otherwise its length is not compatible with VectorSpace
cx = np.array(c[0][:,0], copy=False)
cy = np.array(c[0][:,1], copy=False)
fix = vi.element(cx)
fiy = vi.element(cy)
xp = linspace(0,1,1025)
plot(fx(xp), fy(xp))
plot(fix(xp), fiy(xp))
axis('equal')
Out[17]:
This does not look too good... Let's refine to some more nodes...
In [18]:
Nreps = 11
Nsamples = 100
vi = IteratedVectorSpace(vs, linspace(0,1,Nreps))
x = linspace(0,1,Nsamples)
M = InterpolationMatrix(vi, x)
g = asmatrix(array([fx(x), fy(x)]))
c = lstsq(M, g.T)
cx = np.array(c[0][:,0], copy=False)
cy = np.array(c[0][:,1], copy=False)
fix = vi.element(cx)
fiy = vi.element(cy)
xp = linspace(0,1,1025)
plot(fx(xp), fy(xp))
plot(fix(xp), fiy(xp))
axis('equal')
Out[18]:
In [19]:
q = linspace(0,1,4)
In [20]:
q[2] = -q[2]
In [21]:
q[3] = -q[3]
In [22]:
vs=UniformLagrangeVectorSpace(4)
f = vs.element(q)
fd = vs.element_der(q, 1)
In [25]:
q = [np.array([10,2]), np.array([14,3]),np.array([2,-3])]
print (q)
vs=UniformLagrangeVectorSpace(3)
vs.print_info()
f=vs.element(q)
print (f(4))
In [ ]:
def bisec(f,a,b):
assert(f(a)*f(b) < 0)
toll = 1e-6
err = b - a
n=0
while(err > toll):
c = (b - a) / 2
if(a*c) < 0:
c = b
err = b - a
else:
c = a
err = b - a
n += 1
return [c,n]
a = 0
b = 1
k=bisec(f,a,b)
print (f(a))
In [ ]: