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$

In [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 
from random 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$:

In [2]:
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 eigen-value, 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}$:

$$ \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 $$

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

$$ a y = \lambda m y $$

with matrices $$ 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) $$

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:

In [3]:
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

In [4]:
freedofs = fes.FreeDofs()    
for i in range(len(u.vec)):
    u.vec[i] = random() if freedofs[i] else 0
In [5]:
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(asmall, b=msmall)
    # print (eval / pi**2)
    # print (evec)
    u.vec.data = float(evec[0,0]) * u.vec + float(evec[1,0]) * w
    
Draw (u)

22.751217384109665
2.0881598047601106
2.0021158858300567
2.00021552749329
2.000037586962663
2.0000085868416297
2.000002058016573
2.00000049977882
2.000000122302957
2.000000030115495
2.0000000074658395
2.0000000018764714
2.0000000004924243
2.0000000001488245
2.0000000000633498
2.0000000000420552
2.000000000036742
2.0000000000354152
2.0000000000350853
2.0000000000350027

Simultaneous iteration for several eigenvalues

Declare GridFunction with multiple components to store several eigenvectors.

In [6]:
num = 5
u = GridFunction(fes, multidim=num)

Create list of help vectors:

In [7]:
r = u.vec.CreateVector()
Av = u.vec.CreateVector()
Mv = u.vec.CreateVector()

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

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

Compute num residuals, and solve small eigenvalue problem on $2 num$-dimensional space

In [8]:
lams = num * [1]

asmall = Matrix(2*num, 2*num)
msmall = Matrix(2*num, 2*num)

for i in range(100):
    
    for j in range(num):
        vecs[j].data = u.vecs[j]
        r.data = a.mat * vecs[j] - lams[j] * m.mat * vecs[j]
        # r.data = 1/Norm(r) * r
        r *= 1/Norm(r)
        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[0:num] = ev[0:num]
    print (i, ":", [lam/pi**2 for lam in lams])
    
    for j in range(num):
        r[:] = 0.0
        for k in range(2*num):
            r.data += float(evec[k,j]) * vecs[k]
        u.vecs[j].data = r
    
Draw (u)

0 : [10.592721647004733, 69.832028326716824, 75.723732209472942, 90.073816555536624, 107.98573863773686]
1 : [2.052674041746847, 6.7421542489127768, 10.517758882218944, 14.526785319675637, 17.133918106560674]
2 : [2.0015833169533468, 5.0895887818223056, 5.7211543648400056, 9.995558561124998, 11.323461801775533]
3 : [2.0002482649762627, 5.006455835611006, 5.0789779687837262, 8.5343279528550209, 10.306532711007151]
4 : [2.0000525076691025, 5.0005689734648406, 5.0072647342771504, 8.147891521926045, 10.108803656901301]
5 : [2.0000118908578881, 5.0000744510859247, 5.0007602286015489, 8.0496466302659595, 10.042321062225863]
6 : [2.0000027534673666, 5.0000133612134032, 5.0000994305708728, 8.0190453849181527, 10.016621247546558]
7 : [2.0000006473068996, 5.0000027312886095, 5.0000160176469111, 8.0078602728165915, 10.006548387228728]
8 : [2.0000001535233589, 5.0000005798149623, 5.0000030389299583, 8.0033760282295781, 10.002583797786125]
9 : [2.000000036636385, 5.0000001307917232, 5.000000643744853, 8.0014764928940849, 10.001020896871619]
10 : [2.0000000088004426, 5.0000000331568648, 5.0000001478479357, 8.0006513485445598, 10.000403752197972]
11 : [2.0000000021402342, 5.0000000111235927, 5.0000000378959291, 8.0002885258255088, 10.000159869232109]
12 : [2.0000000005418177, 5.0000000059965126, 5.0000000125272628, 8.0001281204673571, 10.000063386612064]
13 : [2.0000000001572005, 5.0000000047721525, 5.0000000065551298, 8.0000569710711336, 10.000025196845524]
14 : [2.0000000000644897, 5.0000000044561181, 5.0000000051501017, 8.0000253683031382, 10.00001006897358]
15 : [2.000000000042109, 5.0000000043623594, 5.0000000048304489, 8.0000113133902655, 10.000004074256138]
16 : [2.0000000000366995, 5.0000000043347574, 5.0000000047576778, 8.0000050596631294, 10.000001697414231]
17 : [2.0000000000353921, 5.0000000043278527, 5.000000004740599, 8.0000022751867856, 10.000000754771667]
18 : [2.0000000000350751, 5.0000000043261217, 5.0000000047362976, 8.0000010351063011, 10.000000380769581]
19 : [2.0000000000349991, 5.0000000043256732, 5.0000000047353428, 8.0000004825832658, 10.000000232352281]
20 : [2.0000000000349809, 5.0000000043254467, 5.0000000047349591, 8.0000002363675513, 10.000000173435723]
21 : [2.0000000000349751, 5.0000000043255337, 5.0000000047350515, 8.00000012661258, 10.000000150044851]
22 : [2.0000000000349751, 5.0000000043252628, 5.0000000047349626, 8.0000000776820563, 10.000000140755819]
23 : [2.0000000000349742, 5.0000000043254937, 5.0000000047350053, 8.0000000558624187, 10.000000137066214]
24 : [2.0000000000349751, 5.0000000043256083, 5.0000000047355204, 8.0000000461316247, 10.000000135600818]
25 : [2.0000000000349734, 5.0000000043254156, 5.0000000047347735, 8.0000000417910524, 10.000000135018052]
26 : [2.0000000000349738, 5.0000000043250497, 5.0000000047348889, 8.0000000398548146, 10.000000134787244]
27 : [2.0000000000349734, 5.0000000043254786, 5.0000000047348614, 8.0000000389909118, 10.000000134694814]
28 : [2.0000000000349738, 5.0000000043255399, 5.0000000047349102, 8.0000000386054761, 10.000000134658409]
29 : [2.0000000000349751, 5.0000000043252824, 5.0000000047352122, 8.0000000384333543, 10.000000134643781]
30 : [2.0000000000349742, 5.0000000043249253, 5.0000000047348534, 8.0000000383566743, 10.000000134638016]
31 : [2.0000000000349747, 5.0000000043254991, 5.0000000047351465, 8.0000000383224439, 10.00000013463562]
32 : [2.0000000000349738, 5.0000000043252202, 5.0000000047351909, 8.0000000383071228, 10.000000134634966]
33 : [2.0000000000349738, 5.0000000043253872, 5.0000000047349102, 8.0000000383002039, 10.000000134634428]
34 : [2.0000000000349747, 5.0000000043252335, 5.0000000047351225, 8.0000000382971717, 10.000000134634478]
35 : [2.0000000000349742, 5.0000000043252992, 5.0000000047349769, 8.0000000382958358, 10.000000134634419]
36 : [2.0000000000349751, 5.0000000043253543, 5.0000000047347708, 8.0000000382952852, 10.000000134634289]
37 : [2.0000000000349738, 5.0000000043255781, 5.0000000047349147, 8.0000000382950578, 10.000000134633995]
38 : [2.0000000000349729, 5.0000000043253126, 5.0000000047351154, 8.0000000382948908, 10.000000134634293]
39 : [2.0000000000349729, 5.0000000043254218, 5.0000000047348987, 8.000000038294754, 10.000000134634146]
40 : [2.0000000000349738, 5.0000000043254911, 5.0000000047348143, 8.0000000382947292, 10.000000134634289]
41 : [2.0000000000349729, 5.0000000043250976, 5.0000000047346767, 8.000000038294802, 10.000000134634181]
42 : [2.0000000000349742, 5.000000004325253, 5.0000000047348356, 8.0000000382946901, 10.000000134634275]
43 : [2.0000000000349742, 5.0000000043250781, 5.0000000047351545, 8.0000000382945586, 10.000000134634288]
44 : [2.0000000000349742, 5.0000000043251136, 5.0000000047355382, 8.0000000382950898, 10.000000134634535]
45 : [2.0000000000349747, 5.0000000043252131, 5.0000000047348996, 8.0000000382948269, 10.000000134634178]
46 : [2.0000000000349751, 5.0000000043254804, 5.0000000047350612, 8.0000000382947345, 10.000000134634169]
47 : [2.0000000000349747, 5.000000004325452, 5.0000000047350595, 8.0000000382948375, 10.000000134634178]
48 : [2.0000000000349742, 5.0000000043254138, 5.0000000047350079, 8.0000000382948002, 10.000000134634291]
49 : [2.0000000000349738, 5.0000000043256616, 5.0000000047349449, 8.0000000382947984, 10.000000134634337]
50 : [2.0000000000349747, 5.000000004325404, 5.0000000047348676, 8.0000000382947185, 10.000000134634213]
51 : [2.0000000000349751, 5.0000000043254031, 5.0000000047350408, 8.0000000382947825, 10.000000134634325]
52 : [2.0000000000349742, 5.0000000043254627, 5.0000000047350248, 8.0000000382947771, 10.000000134634249]
53 : [2.0000000000349738, 5.0000000043254733, 5.0000000047349191, 8.000000038294818, 10.000000134634302]
54 : [2.0000000000349742, 5.0000000043254147, 5.0000000047350071, 8.0000000382947842, 10.000000134634321]
55 : [2.0000000000349751, 5.0000000043254378, 5.0000000047350159, 8.0000000382947594, 10.00000013463421]
56 : [2.0000000000349751, 5.0000000043254449, 5.0000000047350168, 8.0000000382947505, 10.000000134634263]
57 : [2.0000000000349742, 5.0000000043254369, 5.0000000047350133, 8.0000000382947309, 10.000000134634302]
58 : [2.0000000000349738, 5.0000000043254706, 5.0000000047349644, 8.0000000382947594, 10.000000134634186]
59 : [2.0000000000349751, 5.0000000043253925, 5.0000000047349804, 8.0000000382947007, 10.000000134634249]
60 : [2.0000000000349751, 5.0000000043253987, 5.0000000047349671, 8.0000000382947611, 10.000000134634291]
61 : [2.0000000000349751, 5.0000000043254378, 5.0000000047349822, 8.0000000382947611, 10.000000134634268]
62 : [2.0000000000349742, 5.0000000043253809, 5.0000000047349484, 8.0000000382947078, 10.000000134634215]
63 : [2.0000000000349738, 5.000000004325325, 5.0000000047349422, 8.0000000382948144, 10.000000134634263]
64 : [2.0000000000349747, 5.0000000043254174, 5.0000000047350381, 8.0000000382947842, 10.000000134634252]
65 : [2.0000000000349751, 5.0000000043254191, 5.0000000047349697, 8.0000000382949779, 10.000000134634238]
66 : [2.0000000000349751, 5.000000004325349, 5.0000000047349111, 8.0000000382947398, 10.000000134634252]
67 : [2.0000000000349751, 5.0000000043254591, 5.0000000047349866, 8.0000000382947274, 10.000000134634179]
68 : [2.0000000000349742, 5.0000000043253445, 5.0000000047349173, 8.0000000382947523, 10.000000134634238]
69 : [2.0000000000349738, 5.0000000043254076, 5.0000000047349973, 8.0000000382948002, 10.000000134634332]
70 : [2.0000000000349747, 5.0000000043253809, 5.0000000047349795, 8.0000000382947078, 10.000000134634199]
71 : [2.0000000000349742, 5.0000000043254227, 5.0000000047349449, 8.0000000382947611, 10.000000134634226]
72 : [2.0000000000349747, 5.000000004325396, 5.0000000047350479, 8.0000000382947576, 10.000000134634314]
73 : [2.0000000000349751, 5.000000004325476, 5.0000000047350186, 8.0000000382947327, 10.000000134634254]
74 : [2.0000000000349742, 5.0000000043253818, 5.0000000047349689, 8.0000000382947114, 10.000000134634247]
75 : [2.0000000000349742, 5.0000000043253658, 5.0000000047349964, 8.0000000382947931, 10.000000134634259]
76 : [2.0000000000349751, 5.0000000043254031, 5.0000000047350017, 8.0000000382948144, 10.000000134634313]
77 : [2.0000000000349747, 5.0000000043253348, 5.0000000047349804, 8.0000000382947292, 10.000000134634227]
78 : [2.0000000000349742, 5.0000000043253294, 5.0000000047349848, 8.0000000382947611, 10.000000134634254]
79 : [2.0000000000349751, 5.000000004325476, 5.0000000047350275, 8.0000000382948375, 10.000000134634302]
80 : [2.0000000000349742, 5.0000000043254333, 5.0000000047350026, 8.0000000382948038, 10.000000134634343]
81 : [2.0000000000349747, 5.0000000043253721, 5.0000000047349404, 8.0000000382947931, 10.000000134634275]
82 : [2.0000000000349747, 5.0000000043254502, 5.0000000047350293, 8.0000000382947629, 10.000000134634236]
83 : [2.0000000000349747, 5.0000000043254484, 5.0000000047349982, 8.0000000382947185, 10.000000134634206]
84 : [2.0000000000349742, 5.0000000043254289, 5.0000000047349529, 8.0000000382947167, 10.000000134634254]
85 : [2.0000000000349747, 5.0000000043254236, 5.0000000047350222, 8.0000000382947292, 10.000000134634291]
86 : [2.0000000000349751, 5.0000000043254111, 5.0000000047349937, 8.000000038294738, 10.000000134634252]
87 : [2.0000000000349747, 5.0000000043254289, 5.0000000047349502, 8.00000003829477, 10.000000134634268]
88 : [2.0000000000349742, 5.00000000432542, 5.0000000047349848, 8.0000000382947505, 10.000000134634242]
89 : [2.0000000000349747, 5.0000000043254031, 5.000000004734968, 8.0000000382947274, 10.000000134634275]
90 : [2.0000000000349751, 5.0000000043253205, 5.0000000047349911, 8.0000000382947505, 10.000000134634259]
91 : [2.0000000000349747, 5.0000000043253703, 5.0000000047349005, 8.0000000382947452, 10.000000134634245]
92 : [2.0000000000349751, 5.0000000043253694, 5.0000000047349573, 8.0000000382947611, 10.000000134634238]
93 : [2.0000000000349738, 5.0000000043253925, 5.0000000047349209, 8.0000000382947292, 10.000000134634226]
94 : [2.0000000000349742, 5.0000000043253054, 5.0000000047349795, 8.0000000382947452, 10.000000134634263]
95 : [2.0000000000349738, 5.0000000043253712, 5.000000004734976, 8.0000000382947452, 10.000000134634215]
96 : [2.0000000000349742, 5.0000000043254111, 5.0000000047349644, 8.00000003829485, 10.000000134634313]
97 : [2.0000000000349747, 5.0000000043254218, 5.0000000047351039, 8.0000000382947665, 10.000000134634275]
98 : [2.0000000000349729, 5.0000000043253783, 5.0000000047349982, 8.0000000382947896, 10.000000134634261]
99 : [2.0000000000349747, 5.0000000043254316, 5.0000000047349795, 8.0000000382947789, 10.000000134634258]

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

In [9]: