FreeFEM Documentation on GitHub

stars - forks

Eigen value problems

This section depends on your installation of FreeFEM; you need to have compiled ARPACK. This tool is available in FreeFEM if the word eigenvalue appears in line Load:, like:

1-- FreeFem++ v*.** (date *** *** ** **:**:** CET ****)
2 file : ***.edp
3 Load: lg_fem lg_mesh eigenvalue

This tool is based on arpack++, the object-oriented version of ARPACK eigenvalue package [LEHOUCQ1998].

The function EigenValue computes the generalized eigenvalue of \(A u = \lambda B u\). The Shift-invert method is used by default, with sigma =\(\sigma\) the shift of the method.

The matrix \(OP\) is defined with \(A - \sigma B\).

The return value is the number of converged eigenvalues (can be greater than the number of requested eigenvalues nev=)

1int k = EigenValue(OP, B, nev=Nev, sigma=Sigma);

where the matrix \(OP= A - \sigma B\) with a solver and boundary condition, and the matrix \(B\).

There is also a functional interface:

1int k = EigenValue(n, FOP1, FB, nev=Nev, sigma=Sigma);

where \(n\) is the size of the problem, and the operators are now defined through functions, defining respectively the matrix product of \(OP^{-1}\) and \(B\), as in

1int n = OP1.n;
2func real[int] FOP1(real[int] & u){ real[int] Au = OP^-1*u; return Au; }
3func real[int] FB(real[int] & u){ real[int] Au = B*u; return Au; }

If you want finer control over the method employed in ARPACK, you can specify which mode ARPACK will work with (mode= , see ARPACK documentation [LEHOUCQ1998]). The operators necessary for the chosen mode can be passed through the optional parameters A=, A1=, B=, B1=, (see below).

  • mode=1: Regular mode for solving \(A u = \lambda u\)

    1int k = EigenValue(n, A=FOP, mode=1, nev=Nev);
    

    where the function FOP defines the matrix product of A

  • mode=2: Regular inverse mode for solving \(A u = \lambda B u\)

    1int k = EigenValue(n, A=FOP, B=FB, B1=FB1, mode=2, nev=Nev);
    

    where the functions FOP, FB and FB1 define respectively the matrix product of \(A\), \(B\) and \(B^{-1}\)

  • mode=3: Shift-invert mode for solving \(A u = \lambda B u\)

    1int k = EigenValue(n, A1=FOP1, B=FB, mode=3, sigma=Sigma, nev=Nev);
    

    where the functions FOP1 and FB define respectively the matrix product of \(OP^{-1} = (A - \sigma B)^{-1}\) and \(B\)

You can also specify which subset of eigenvalues you want to compute (which=). The default value is which="LM", for eigenvalues with largest magnitude. "SM" is for smallest magnitude, "LA" for largest algebraic value, "SA" for smallest algebraic value, and "BE" for both ends of the spectrum.

Remark: For complex problems, you need to use the keyword complexEigenValue instead of EigenValue when passing operators through functions.

Note

Boundary condition and Eigenvalue Problems

The locking (Dirichlet) boundary condition is make with exact penalization so we put 1e30=tgv on the diagonal term of the locked degree of freedom (see Finite element chapter). So take Dirichlet boundary condition just on \(A\) and not on \(B\) because we solve \(w=OP^{-1}*B*v\).

If you put locking (Dirichlet) boundary condition on \(B\) matrix (with key work on) you get small spurious modes \((10^{-30})\), due to boundary condition, but if you forget the locking boundary condition on \(B\) matrix (no keywork on) you get huge spurious \((10^{30})\) modes associated to these boundary conditons. We compute only small mode, so we get the good one in this case.

  • sym= The problem is symmetric (all the eigen value are real)

  • nev= The number desired eigenvalues (nev) close to the shift.

  • value= The array to store the real part of the eigenvalues

  • ivalue= The array to store the imaginary part of the eigenvalues

  • vector= The FE function array to store the eigenvectors

  • rawvector= An array of type real[int,int] to store eigenvectors by column.

    For real non symmetric problems, complex eigenvectors are given as two consecutive vectors, so if eigenvalue \(k\) and \(k+1\) are complex conjugate eigenvalues, the \(k\)th vector will contain the real part and the \(k+1\)th vector the imaginary part of the corresponding complex conjugate eigenvectors.

  • tol= The relative accuracy to which eigenvalues are to be determined;

  • sigma= The shift value;

  • maxit= The maximum number of iterations allowed;

  • ncv= The number of Arnoldi vectors generated at each iteration of ARPACK;

  • mode= The computational mode used by ARPACK (see above);

  • which= The requested subset of eigenvalues (see above).

Tip

Laplace eigenvalue

In the first example, we compute the eigenvalues and the eigenvectors of the Dirichlet problem on square \(\Omega=]0,\pi[^2\).

The problem is to find: \(\lambda\), and \(\nabla u_{\lambda}\) in \(\mathbb{R}{\times} H^1_0(\Omega)\)

\[\int_\Omega \nabla u_{\lambda} \nabla v = \lambda \int_\Omega u v \quad \forall v \in H^1_0(\Omega)\]

The exact eigenvalues are \(\lambda_{n,m} =(n^2+m^2), (n,m)\in {\mathbb{N}_*}^2\) with the associated eigenvectors are \(u_{{m,n}}=\sin(nx)*\sin(my)\).

We use the generalized inverse shift mode of the arpack++ library, to find 20 eigenvalues and eigenvectors close to the shift value \(\sigma=20\).

 1// Parameters
 2verbosity=0;
 3real sigma = 20; //value of the shift
 4int nev = 20; //number of computed eigen value close to sigma
 5
 6// Mesh
 7mesh Th = square(20, 20, [pi*x, pi*y]);
 8
 9// Fespace
10fespace Vh(Th, P2);
11Vh u1, u2;
12
13// Problem
14// OP = A - sigma B ; // the shifted matrix
15varf op (u1, u2)
16    = int2d(Th)(
17          dx(u1)*dx(u2)
18        + dy(u1)*dy(u2)
19        - sigma* u1*u2
20    )
21    + on(1, 2, 3, 4, u1=0)
22    ;
23
24varf b ([u1], [u2]) = int2d(Th)(u1*u2); //no boundary condition
25
26matrix OP = op(Vh, Vh, solver=Crout, factorize=1); //crout solver because the matrix in not positive
27matrix B = b(Vh, Vh, solver=CG, eps=1e-20);
28
29// important remark:
30// the boundary condition is make with exact penalization:
31// we put 1e30=tgv on the diagonal term of the lock degree of freedom.
32// So take Dirichlet boundary condition just on $a$ variational form
33// and not on $b$ variational form.
34// because we solve $ w=OP^-1*B*v $
35
36// Solve
37real[int] ev(nev); //to store the nev eigenvalue
38Vh[int] eV(nev); //to store the nev eigenvector
39
40int k = EigenValue(OP, B, sym=true, sigma=sigma, value=ev, vector=eV,
41    tol=1e-10, maxit=0, ncv=0);
42
43// Display & Plot
44for (int i = 0; i < k; i++){
45    u1 = eV[i];
46    real gg = int2d(Th)(dx(u1)*dx(u1) + dy(u1)*dy(u1));
47    real mm = int2d(Th)(u1*u1) ;
48    cout << "lambda[" << i << "] = " << ev[i] << ", err= " << int2d(Th)(dx(u1)*dx(u1) + dy(u1)*dy(u1) - (ev[i])*u1*u1) << endl;
49    plot(eV[i], cmm="Eigen Vector "+i+" value ="+ev[i], wait=true, value=true);
50}

The output of this example is:

 1lambda[0] = 5.0002, err= -1.46519e-11
 2lambda[1] = 8.00074, err= -4.05158e-11
 3lambda[2] = 10.0011, err= 2.84925e-12
 4lambda[3] = 10.0011, err= -7.25456e-12
 5lambda[4] = 13.002, err= -1.74257e-10
 6lambda[5] = 13.0039, err= 1.22554e-11
 7lambda[6] = 17.0046, err= -1.06274e-11
 8lambda[7] = 17.0048, err= 1.03883e-10
 9lambda[8] = 18.0083, err= -4.05497e-11
10lambda[9] = 20.0096, err= -2.21678e-13
11lambda[10] = 20.0096, err= -4.16212e-14
12lambda[11] = 25.014, err= -7.42931e-10
13lambda[12] = 25.0283, err= 6.77444e-10
14lambda[13] = 26.0159, err= 3.19864e-11
15lambda[14] = 26.0159, err= -4.9652e-12
16lambda[15] = 29.0258, err= -9.99573e-11
17lambda[16] = 29.0273, err= 1.38242e-10
18lambda[17] = 32.0449, err= 1.2522e-10
19lambda[18] = 34.049, err= 3.40213e-11
20lambda[19] = 34.0492, err= 2.41751e-10
EigenValueProblems1

Fig. 159 Isovalue of 11th eigenvector \(u_{4,3}-u_{3,4}\)

EigenValueProblems2

Fig. 160 Isovalue of 12th eigenvector \(u_{4,3}+u_{3,4}\)

Table of content