Geometric Multigrid

In [1]:

    

%matplotlib notebook
import numpy

m = 101
A = numpy.eye(m) - .5*numpy.eye(m,k=1) - .5*numpy.eye(m,k=-1)

I = numpy.eye(m)
L, X  = numpy.linalg.eigh(A)

from matplotlib import pyplot
pyplot.style.use('ggplot')
pyplot.figure()
pyplot.plot(L, '.')
pyplot.show()

print(L[:4], L[-4:])
print('cond', L[-1]/L[0])


    

    















    












    




[ 0.00047428  0.00189667  0.00426582  0.00757949] [ 1.99242051  1.99573418  1.99810333  1.99952572]
cond 4215.9157276

In [2]:

    

pyplot.figure()
for i in (0,1,2):
    pyplot.plot(X[:,i], label='$\lambda = %f$' % L[i])
pyplot.legend(loc='upper right')
pyplot.show()


    

In [3]:

    

pyplot.figure()
for i in (10,90):
    pyplot.plot(X[:,i], '-o', label='$\lambda = %f$' % L[i])
pyplot.legend(loc='upper right')
pyplot.show()


    

In [4]:

    

P = numpy.eye(m) + .5*numpy.eye(m,k=1) + .5*numpy.eye(m,k=-1)
P = P[:,::2]
#P = numpy.eye(m) + 2./3*(numpy.eye(m,k=1)+numpy.eye(m,k=-1)) + \
#    1./3*(numpy.eye(m,k=2)+numpy.eye(m,k=-2))
#P = P[:,::3]

pyplot.figure()
pyplot.plot(P[:10,2:4], '-*')

P.shape


    

    















    












    
Out[4]:





(101, 51)

In [5]:

    

pyplot.figure()
pyplot.plot(1./2*P.dot(P.T.dot(X[:,5])), label='coarse')
pyplot.plot(X[:,5], label='fine')
pyplot.legend(loc='upper right')


    

    















    












    
Out[5]:





<matplotlib.legend.Legend at 0x7f08411da198>

In [6]:

    

Ac = P.T.dot(A.dot(P))
print(Ac.shape)

Lc, Xc = numpy.linalg.eigh(Ac)
pyplot.figure()
pyplot.plot(Lc, '*')


    

    




(51, 51)






    















    












    
Out[6]:





[<matplotlib.lines.Line2D at 0x7f084114c2e8>]

In [7]:

    

pyplot.figure()
pyplot.spy(Ac)
Ac[:5,:5]


    

    















    












    
Out[7]:





array([[ 0.75, -0.25,  0.  ,  0.  ,  0.  ],
       [-0.25,  0.5 , -0.25,  0.  ,  0.  ],
       [ 0.  , -0.25,  0.5 , -0.25,  0.  ],
       [ 0.  ,  0.  , -0.25,  0.5 , -0.25],
       [ 0.  ,  0.  ,  0.  , -0.25,  0.5 ]])

In [8]:

    

# Consider the A-orthogonal projector into the range of P
Sc = P.dot(numpy.linalg.inv(Ac)).dot(P.T).dot(A)
Ls, Xs = numpy.linalg.eig(Sc)
print(max(abs(Ls.imag)))
Ls = Ls.real
idx = Ls.argsort()
Ls = Ls[idx]; Xs = Xs[:,idx]

pyplot.figure()
pyplot.plot(Ls, '.')


    

    




7.18533924784e-15






    















    












    
Out[8]:





[<matplotlib.lines.Line2D at 0x7f08410a2710>]

In [9]:

    

# Iteration matrix for a V(1,1) cycle
Lt, Xt = numpy.linalg.eig((I - .67*A).dot(I - Sc).dot(I - .67*A))
pyplot.figure()
pyplot.plot(Lt.real, Lt.imag, '.')


    

    















    












    
Out[9]:





[<matplotlib.lines.Line2D at 0x7f0841052240>]

In [ ]:

    




    

In [ ]: