This page was generated from jupyter-files/unit-2.2-eigenvalues/pinvit.ipynb.

2.2 Eigenvalue problems

We solve the generalized eigenvalue problem

\[A u = \lambda M u,\]

where \(A\) comes from \(\int \nabla u \nabla v\), and \(M\) from \(\int u v\), on the space \(H_0^1\)

import netgen.gui
%gui tk
from netgen.geom2d import unit_square
from ngsolve import *
from math import pi
import scipy.linalg
from scipy import random

mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))

We setup a stiffness matrix \(A\) and a mass matrix \(M\), and declare a preconditioner for \(A\):

fes = H1(mesh, order=4, dirichlet=[1,2,3,4])
u = fes.TrialFunction()
v = fes.TestFunction()

a = BilinearForm(fes)
a += SymbolicBFI(grad(u)*grad(v))
pre = Preconditioner(a, "multigrid")

m = BilinearForm(fes)
m += SymbolicBFI(u*v)

a.Assemble()
m.Assemble()

u = GridFunction(fes)

The inverse iteration is

\[u_{n+1} = A^{-1} M u_n,\]

where the Rayleigh quotient

\[\rho_n = \frac{\left \langle A u_n, u_n\right \rangle}{\left \langle M u_n, u_n\right \rangle}\]

converges to the smallest eigenvalue, with rate of convergence \(\lambda_1 / \lambda_2\).

The preconditioned inverse iteration (PINVIT), see [Knyazef+Neymeyr], replaces \(A^{-1}\) by an approximate inverse \(C^{-1}\):

\[\begin{split}\rho_n = \frac{\left \langle A u_n, u_n\right \rangle}{\left \langle M u_n, u_n\right \rangle} \\ w_n = C^{-1} (A u_n - \rho M u_n) \\ u_{n+1} = u_n + \alpha w_n\end{split}\]

The optimal step-size \(\alpha\) is found by minimizing the Rayleigh-quotient on a two-dimensional space:

\[u_{n+1} = \operatorname{arg} \min_{v \in \{ u_n, w_n\}} \frac{\left \langle A v, v\right \rangle}{\left \langle M v, v\right \rangle}\]

This minimization problem can be solved by a small eigenvalue problem

\[ \begin{align}\begin{aligned} a y = \lambda m y\\with matrices\end{aligned}\end{align} \]

\[\begin{split}a = \left( \begin{array}{cc} \left \langle A u_n, u_n \right \rangle & \left \langle A u_n, w_n \right \rangle \\ \left \langle A w_n, u_n \right \rangle & \left \langle A w_n, w_n \right \rangle \end{array} \right) \quad m = \left( \begin{array}{cc} \left \langle M u_n, u_n \right \rangle & \left \langle M u_n, w_n \right \rangle \\ \left \langle M w_n, u_n \right \rangle & \left \langle M w_n, w_n \right \rangle \end{array} \right)\end{split}\]

Then, the new iterate is

\[u_{n+1} = y_1 u_n + y_2 w_n\]

where \(y\) is the eigenvector corresponding to the smaller eigenvalue.

Implementation in NGSolve

Create some help vectors:

r = u.vec.CreateVector()
w = u.vec.CreateVector()
Mu = u.vec.CreateVector()
Au = u.vec.CreateVector()
Mw = u.vec.CreateVector()
Aw = u.vec.CreateVector()

random initial condition, zero on the boundary

r.FV().NumPy()[:] = random.rand(fes.ndof)
u.vec.data = Projector(fes.FreeDofs(), True) * r

# from random import random
# freedofs = fes.FreeDofs()
# for i in range(len(u.vec)):
#    u.vec[i] = random() if freedofs[i] else 0

for i in range(20):
    Au.data = a.mat * u.vec
    Mu.data = m.mat * u.vec
    auu = InnerProduct(Au, u.vec)
    muu = InnerProduct(Mu, u.vec)
    # Rayleigh quotiont
    lam = auu/muu
    print (lam / (pi**2))
    # residual
    r.data = Au - lam * Mu
    w.data = pre.mat * r.data
    w.data = 1/Norm(w) * w
    Aw.data = a.mat * w
    Mw.data = m.mat * w

    # setup and solve 2x2 small eigenvalue problem
    asmall = Matrix(2,2)
    asmall[0,0] = auu
    asmall[0,1] = asmall[1,0] = InnerProduct(Au, w)
    asmall[1,1] = InnerProduct(Aw, w)
    msmall = Matrix(2,2)
    msmall[0,0] = muu
    msmall[0,1] = msmall[1,0] = InnerProduct(Mu, w)
    msmall[1,1] = InnerProduct(Mw, w)
    # print ("asmall =", asmall, ", msmall = ", msmall)


    eval,evec = scipy.linalg.eigh(a=asmall, b=msmall)
    # print (eval, evec)
    u.vec.data = float(evec[0,0]) * u.vec + float(evec[1,0]) * w

Draw (u)

26.310136660400275
2.0782284240636826
2.0016788490372455
2.0002002507016474
2.0000413285947802
2.0000105080928825
2.0000027530788227
2.0000007313313577
2.0000001959000486
2.0000000528126267
2.0000000143203
2.000000003918371
2.0000000010972063
2.000000000329983
2.0000000001208567
2.0000000000637526
2.000000000048131
2.000000000043857
2.0000000000426827
2.0000000000423617

Simultaneous iteration for several eigenvalues

Declare GridFunction with multiple components to store several eigenvectors.

num = 5
u = GridFunction(fes, multidim=num)

Create list of help vectors, random initialization of u vectors:

r = u.vec.CreateVector()
Av = u.vec.CreateVector()
Mv = u.vec.CreateVector()

vecs = []
for i in range(2*num):
    vecs.append (u.vec.CreateVector())

for v in u.vecs:
    r.FV().NumPy()[:] = random.rand(fes.ndof)
    v.data = Projector(fes.FreeDofs(), True) * r

Compute num residuals, and solve small eigenvalue problem on \(2 num\)-dimensional space

asmall = Matrix(2*num, 2*num)
msmall = Matrix(2*num, 2*num)
lams = num * [1]

for i in range(20):

    for j in range(num):
        vecs[j].data = u.vecs[j]
        r.data = a.mat * vecs[j] - lams[j] * m.mat * vecs[j]
        vecs[num+j].data = pre.mat * r

    for j in range(2*num):
        Av.data = a.mat * vecs[j]
        Mv.data = m.mat * vecs[j]
        for k in range(2*num):
            asmall[j,k] = InnerProduct(Av, vecs[k])
            msmall[j,k] = InnerProduct(Mv, vecs[k])

    ev,evec = scipy.linalg.eigh(a=asmall, b=msmall)
    lams[:] = ev[0:num]
    print (i, ":", [lam/pi**2 for lam in lams])

    for j in range(num):
        u.vecs[j][:] = 0.0
        for k in range(2*num):
            u.vecs[j].data += float(evec[k,j]) * vecs[k]

Draw (u)

0 : [11.6130455768801, 61.066110947635224, 82.49711365113266, 90.41384577806575, 96.98583176269047]
1 : [2.037296014502914, 5.994129376899316, 9.212270101228079, 16.752145334653957, 19.356324864358374]
2 : [2.0014100615066996, 5.092284538737162, 5.238911096053245, 10.123988712254931, 10.879020322598242]
3 : [2.0002554163633968, 5.012268562196962, 5.016722138616416, 9.213915664378316, 10.256909163977726]
4 : [2.000057169527135, 5.001025982274766, 5.0025582296087965, 8.719831456835285, 10.0918178567872]
5 : [2.000013593722611, 5.000120805946273, 5.000384594365776, 8.377992268553804, 10.035379830359405]
6 : [2.000003314919647, 5.000020247177596, 5.000055876284345, 8.184287429684003, 10.01410422992045]
7 : [2.000000821205276, 5.000004016939583, 5.00000849353781, 8.08609184849783, 10.005671478126507]
8 : [2.0000002058391035, 5.0000008179428255, 5.000001478437589, 8.039447631058222, 10.002288011579898]
9 : [2.0000000520461136, 5.000000168593007, 5.000000304198153, 8.017928866313158, 10.000924591937792]
10 : [2.000000013279318, 5.000000039218533, 5.000000071766628, 8.008121745129989, 10.000374075767763]
11 : [2.0000000034302508, 5.000000012511654, 5.000000020495333, 8.003675600991937, 10.000151525029434]
12 : [2.0000000009138295, 5.000000006632997, 5.0000000086298915, 8.00166315351873, 10.00006147486362]
13 : [2.000000000267371, 5.000000005221903, 5.000000005864453, 8.00075273753271, 10.000025013254238]
14 : [2.0000000001006097, 5.000000004822398, 5.000000005260717, 8.00034080330491, 10.000010241080249]
15 : [2.000000000057419, 5.00000000470997, 5.000000005129163, 8.000154364958556, 10.000004253440505]
16 : [2.0000000000461995, 5.000000004681857, 5.000000005098049, 8.00006995469134, 10.000001825435998]
17 : [2.000000000043275, 5.000000004674421, 5.000000005090545, 8.000031724949128, 10.000000840538995]
18 : [2.0000000000425104, 5.000000004672643, 5.000000005088353, 8.000014405452335, 10.000000440900093]
19 : [2.00000000004231, 5.000000004672329, 5.000000005087833, 8.000006556850161, 10.000000278698185]

The multidim-component select in the Visualization dialog box allows to switch between the components of the multidim-GridFunction.