Skip to main content
`Ctrl`+`K`

# ex13.py: Nonlinear eigenproblem with contour integral

This example solves a nonlinear eigenvalue problem with the `NEP` module configured to apply a contour integral method.

The problem arises from the PDE

\[u_{xx}(x) + n_c^2 \lambda^2 u(x) + g(\lambda) D_0 \lambda^2 u(x) = 0\]

where

\[\begin{split}g(\lambda) &= g_t/(\lambda-k_a + i g_t), \\ k_a &= 8.0, \\ g_t &= 0.5, \\ D_0 &= 0.5, \\ n_c &= 1.2,\end{split}\]

and the boundary conditions are

\[\begin{split}u(0) &= 0 \\ u_x(1) &= i \lambda u(1).\end{split}\]

For the discretization, \(n\) grid points are used, \(x_1=0,\dots,x_n=1\), with step size \(h = 1/(n-1)\). Hence,

\[\begin{split}u_{xx}(x_i) &= \frac{1}{h^2} (u_{i-1} - 2 u_i + u_{i+1}) \\ &= \frac{1}{h^2} [1, -2, 1] [u_{i-1}, u_i, u_{i+1}]^T.\end{split}\]

The boundary condition at \(x=0\) is \(u_1=0\), and at \(x=1\):

\[\begin{split}u'(1) &= 1/2 ( u'(1+h/2) + u'(1-h/2) ) \\ &= 1/2 ( (u_{n+1}-u_n)/h + (u_n-u_{n-1})/h ) \\ &= 1/(2h) (u_{n+1} - u_{n-1}) = i\lambda u_n,\end{split}\]

and therefore

\[u_{n+1} = 2i h \lambda u_n + u_{n-1}.\]

The Laplace term for \(u_n\) is

\[\begin{split}&= \frac{1}{h^2} (u_{n-1} - 2u_n + u_{n+1}) \\ &= \frac{1}{h^2} (u_{n-1} - 2u_n + 2ih\lambda u_n + u_{n-1}) \\ &= \frac{1}{h^2} (2 u_{n-1} + (2ih\lambda - 2) u_n) \\ &= \frac{1}{h^2} [2, 2ih\lambda -2] [u_{n-1}, u_n]^T.\end{split}\]

The above discretization allows us to write the nonlinear PDE in the following split-operator form

\[\{A + \lambda^2 n_c^2 I_d + g(\lambda) \lambda^2 D_0 I_d + 2i\lambda/h D\} u = 0\]

so \(f_1 = 1\), \(f_2 = n_c^2 \lambda^2\), \(f_3 = g(\lambda) \lambda^2D_0\), \(f_4 = 2i\lambda/h\), with coefficient matrices

\[\begin{split}A = \begin{bmatrix} 1 & 0 & 0 & \dots & 0 \\ 0 & * & * & \dots & 0 \\ 0 & * & * & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & * \end{bmatrix}, I_d = \begin{bmatrix} 0 & 0 & 0 & \dots & 0 \\ 0 & 1 & 0 & \dots & 0 \\ 0 & 0 & 1 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & 0 \end{bmatrix}, D = \begin{bmatrix} 0 & 0 & 0 & \dots & 0 \\ 0 & 0 & 0 & \dots & 0 \\ 0 & 0 & 0 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & 1 \end{bmatrix}.\end{split}\]

Contributed by: Thomas Hisch.

The full source code for this demo can be downloaded here.

Initialization by importing slepc4py, petsc4py, numpy and scipy.

    import sys

    import slepc4py

    slepc4py.init(sys.argv)  # isort:skip

    import numpy as np

    try:
        import scipy
        import scipy.optimize
    except ImportError:
        scipy = None

    from petsc4py import PETSc
    from slepc4py import SLEPc

    Print = PETSc.Sys.Print

Check that the selected PETSc/SLEPc are built with complex scalars.

    if not np.issubdtype(PETSc.ScalarType, np.complexfloating):
        Print("Demo should only be executed with complex PETSc scalars")
        exit(0)

Long function that defines the nonlinear eigenproblem and solves it.

    def solve(n):
        L = 1.0
        h = L / (n - 1)
        nc = 1.2
        ka = 10.0
        gt = 4.0
        D0 = 0.5

        A = PETSc.Mat().create()
        A.setSizes([n, n])
        A.setFromOptions()
        A.setOption(PETSc.Mat.Option.HERMITIAN, False)

        rstart, rend = A.getOwnershipRange()
        d0, d1, d2 = (
            1 / h**2,
            -2 / h**2,
            1 / h**2,
        )
        Print(f"dterms={(d0, d1, d2)}")

        if rstart == 0:
            # dirichlet boundary condition at the left lead
            A[0, 0] = 1.0
            A[0, 1] = 0.0
            A[1, 0] = 0.0

            A[1, 1] = d1
            A[1, 2] = d2
            rstart += 2
        if rend == n:
            # at x=1.0 neumann boundary condition (not handled here but in a
            # different matrix (D))
            A[n - 1, n - 2] = 2.0 / h**2
            A[n - 1, n - 1] = (-2) / h**2  # + 2j*k*h / h**2 (neumann)
            rend -= 1

        for i in range(rstart, rend):
            A[i, i - 1 : i + 2] = [d0, d1, d2]

        A.assemble()

        Id = PETSc.Mat().create()
        Id.setSizes([n, n])
        Id.setFromOptions()
        Id.setOption(PETSc.Mat.Option.HERMITIAN, True)
        rstart, rend = Id.getOwnershipRange()
        if rstart == 0:
            # due to dirichlet BC
            rstart += 1
        for i in range(rstart, rend):
            Id[i, i] = 1.0
        Id.assemble()

        D = PETSc.Mat().create()
        D.setSizes([n, n])
        D.setFromOptions()
        D.setOption(PETSc.Mat.Option.HERMITIAN, True)
        _, rend = D.getOwnershipRange()
        if rend == n:
            D[n - 1, n - 1] = 1
        D.assemble()

        Print(f"DOF: {A.getInfo()['nz_used']}, MEM: {A.getInfo()['memory']}")

        f1 = SLEPc.FN().create()
        f1.setType(SLEPc.FN.Type.RATIONAL)
        f1.setRationalNumerator([1.0])
        f2 = SLEPc.FN().create()
        f2.setType(SLEPc.FN.Type.RATIONAL)
        f2.setRationalNumerator([nc**2, 0.0, 0.0])
        f3 = SLEPc.FN().create()
        f3.setType(SLEPc.FN.Type.RATIONAL)
        f3.setRationalNumerator([D0 * gt, 0.0, 0.0])
        f3.setRationalDenominator([1.0, -ka + 1j * gt])
        f4 = SLEPc.FN().create()
        f4.setType(SLEPc.FN.Type.RATIONAL)
        f4.setRationalNumerator([2j / h, 0])

        # Setup the solver
        nep = SLEPc.NEP().create()

        nep.setSplitOperator(
            [A, Id, Id, D],
            [f1, f2, f3, f4],
            PETSc.Mat.Structure.SUBSET,
        )

        # Customize options
        nep.setTolerances(tol=1e-7)

        nep.setDimensions(nev=24)
        nep.setType(SLEPc.NEP.Type.CISS)

        # the rg params are chosen s.t. the singularity at k = ka - 1j*gt is
        # outside of the contour.
        radius = 3 * gt
        vscale = 0.5 * gt / radius
        rg_params = (ka, 3 * gt, vscale)
        R = nep.getRG()
        R.setType(SLEPc.RG.Type.ELLIPSE)
        Print(f"RG params: {rg_params}")
        R.setEllipseParameters(*rg_params)

        nep.setFromOptions()

        # Solve the problem
        nep.solve()

        its = nep.getIterationNumber()
        Print("Number of iterations of the method: %i" % its)
        sol_type = nep.getType()
        Print("Solution method: %s" % sol_type)
        nev, ncv, mpd = nep.getDimensions()
        Print("")
        Print("Subspace dimension: %i" % ncv)
        tol, maxit = nep.getTolerances()
        Print("Stopping condition: tol=%.4g" % tol)
        Print("")

        nconv = nep.getConverged()
        Print("Number of converged eigenpairs %d" % nconv)

        x = A.createVecs("right")

        evals = []
        modes = []
        if nconv > 0:
            Print()
            Print("        lam            ||T(lam)x||   |lam-lam_exact|/|lam_exact| ")
            Print("--------------------- ------------- -----------------------------")
            for i in range(nconv):
                lam = nep.getEigenpair(i, x)
                error = nep.computeError(i)

                def eigenvalue_error_term(k):
                    gkmu = gt / (k - ka + 1j * gt)
                    nceff = np.sqrt(nc**2 + gkmu * D0)
                    return -1j / np.tan(nceff * k * L) - 1 / nceff

                # compute the expected_eigenvalue

                # we assume that the numerically calculated eigenvalue is close to
                # the exact one, which we can determine using a Newton-Raphson
                # method.
                if scipy:
                    expected_lam = scipy.optimize.newton(
                        eigenvalue_error_term, np.complex128(lam), rtol=1e-11
                    )
                    rel_err = abs(lam - expected_lam) / abs(expected_lam)
                    rel_err = "%6g" % rel_err
                else:
                    rel_err = "scipy not installed"

                Print(" %9f%+9f j %12g   %s" % (lam.real, lam.imag, error, rel_err))

                evals.append(lam)
                modes.append(x.getArray().copy())
            Print()

        return np.asarray(evals), rg_params, ka, gt

The main function reads the problem size `n` from the command line, solves the problem with the above function, and then plots the computed eigenvalues.

    def main():
        opts = PETSc.Options()
        n = opts.getInt("n", 256)
        Print(f"n={n}")

        evals, rg_params, ka, gt = solve(n)

        if not opts.getBool("ploteigs", True) or PETSc.COMM_WORLD.getRank():
            return

        try:
            import matplotlib.pyplot as plt
            from matplotlib.patches import Ellipse
        except ImportError:
            print("plot is not shown, because matplotlib is not installed")
        else:
            fig, ax = plt.subplots()
            ax.plot(evals.real, evals.imag, "x")

            height = 2 * rg_params[1] * rg_params[2]
            ellipse = Ellipse(
                xy=(rg_params[0], 0.0),
                width=rg_params[1] * 2,
                height=height,
                edgecolor="r",
                fc="None",
                lw=2,
            )
            ax.add_patch(ellipse)

            ax.grid()
            ax.legend()
            plt.show()

    if __name__ == "__main__":
        main()