Skip to main content
`Ctrl`+`K`

# ex7.py: Nonlinear eigenproblem with callback functions

This example solves a nonlinear eigenvalue problem arising from the the discretization of a PDE on a one-dimensional domain with finite differences. The nonlinearity comes from the boundary conditions. The PDE is

\[-u'' = \lambda u\]

defined on the interval [0,1] and subject to the boundary conditions

\[u(0)=0, u'(1)=u(1)\lambda\frac{\kappa}{\kappa-\lambda},\]

where \(\lambda\) is the eigenvalue and \(\kappa\) is a parameter.

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

Initialization is similar to previous examples.

    import sys, slepc4py
    slepc4py.init(sys.argv)

    from petsc4py import PETSc
    from slepc4py import SLEPc
    from numpy import sqrt, sin

    Print = PETSc.Sys.Print

When implementing a nonlinear eigenproblem with callback functions we must provide code that builds the function matrix \(T(\lambda)\) for a given \(\lambda\) and optionally the Jacobian matrix \(T'(\lambda)\), i.e., the derivative with respect to the eigenvalue.

In slepc4py the callbacks are integrated in a class. In this example, apart from the constructor, we have three methods:

- `formFunction` to fill the function matrix `F`. Note that `F` is received as an argument and we just need to fill its entries using the value of the parameter `mu`. Matrix `B` is used to build the preconditioner, and is usually equal to `F`.

- `formJacobian` to fill the Jacobian matrix `J`. Some eigensolvers do not need this, but it is recommended to implement it.

- `checkSolution` is just a convenience method to check that a given solution satisfies the PDE.

    class MyPDE(object):

        def __init__(self, kappa, h):
            self.kappa = kappa
            self.h     = h

        def formFunction(self, nep, mu, F, B):
            n, m = F.getSize()
            Istart, Iend = F.getOwnershipRange()
            i1 = Istart
            if Istart==0: i1 = i1 + 1
            i2 = Iend
            if Iend==n: i2 = i2 - 1
            h = self.h
            c = self.kappa/(mu-self.kappa)
            d = n

            # Interior grid points
            for i in range(i1,i2):
                val = -d-mu*h/6.0
                F[i,i-1] = val
                F[i,i]   = 2.0*(d-mu*h/3.0)
                F[i,i+1] = val

            # Boundary points
            if Istart==0:
                F[0,0] = 2.0*(d-mu*h/3.0)
                F[0,1] = -d-mu*h/6.0
            if Iend==n:
                F[n-1,n-2] = -d-mu*h/6.0
                F[n-1,n-1] = d-mu*h/3.0+mu*c

            F.assemble()
            if B != F: B.assemble()
            return PETSc.Mat.Structure.SAME_NONZERO_PATTERN

        def formJacobian(self, nep, mu, J):
            n, m = J.getSize()
            Istart, Iend = J.getOwnershipRange()
            i1 = Istart
            if Istart==0: i1 = i1 + 1
            i2 = Iend
            if Iend==n: i2 = i2 - 1
            h = self.h
            c = self.kappa/(mu-self.kappa)

            # Interior grid points
            for i in range(i1,i2):
                J[i,i-1] = -h/6.0
                J[i,i]   = -2.0*h/3.0
                J[i,i+1] = -h/6.0

            # Boundary points
            if Istart==0:
                J[0,0] = -2.0*h/3.0
                J[0,1] = -h/6.0
            if Iend==n:
                J[n-1,n-2] = -h/6.0
                J[n-1,n-1] = -h/3.0-c*c

            J.assemble()
            return PETSc.Mat.Structure.SAME_NONZERO_PATTERN

        def checkSolution(self, mu, y):
            nu = sqrt(mu)
            u = y.duplicate()
            n = u.getSize()
            Istart, Iend = J.getOwnershipRange()
            h = self.h
            for i in range(Istart,Iend):
                x = (i+1)*h
                u[i] = sin(nu*x);
            u.assemble()
            u.normalize()
            u.axpy(-1.0,y)
            return u.norm()

We use an auxiliary function `FixSign` to force the computed eigenfunction to be real and positive, since some eigensolvers may return the eigenvector multiplied by a complex number of modulus one.

    def FixSign(x):
        comm = x.getComm()
        rank = comm.getRank()
        n = 1 if rank == 0 else 0
        aux = PETSc.Vec().createMPI((n, PETSc.DECIDE), comm=comm)
        if rank == 0: aux[0] = x[0]
        aux.assemble()
        x0 = aux.sum()
        sign = x0/abs(x0)
        x.scale(1.0/sign)

The main program processes two command-line options, `n` (size of the grid) and `kappa` (the parameter of the PDE), then creates an object of the class we have defined previously.

    opts = PETSc.Options()
    n = opts.getInt('n', 128)
    kappa = opts.getReal('kappa', 1.0)
    pde = MyPDE(kappa, 1.0/n)

In order to set up the solver we have to pass the two callback functions (methods of the class) together with the matrix objects that will be used every time these methods are called. In this simple example we can do a preallocation of the matrices, although this is not necessary.

    nep = SLEPc.NEP().create()
    F = PETSc.Mat().create()
    F.setSizes([n, n])
    F.setType('aij')
    F.setPreallocationNNZ(3)
    nep.setFunction(pde.formFunction, F)

    J = PETSc.Mat().create()
    J.setSizes([n, n])
    J.setType('aij')
    J.setPreallocationNNZ(3)
    nep.setJacobian(pde.formJacobian, J)

After setting some options, we can solve the problem. Here we also illustrate how to pass an initial guess to the solver.

    nep.setTolerances(tol=1e-9)
    nep.setDimensions(1)
    nep.setFromOptions()

    x = F.createVecs('right')
    x.set(1.0)
    nep.setInitialSpace(x)
    nep.solve()

Once the solver has finished, we print some information together with the computed solution. For each computed eigenpair, we print the residual norm and also the error estimated with the class method `checkSolution`.

    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 )

    if nconv > 0:
      Print()
      Print("        k              ||T(k)x||          error ")
      Print("----------------- ------------------ ------------------")
      for i in range(nconv):
        k = nep.getEigenpair(i, x)
        FixSign(x)
        res = nep.computeError(i)
        error = pde.checkSolution(k.real, x)
        if k.imag != 0.0:
          Print( " %9f%+9f j %12g     %12g" % (k.real, k.imag, res, error) )
        else:
          Print( " %12f       %12g     %12g" % (k.real, res, error) )
      Print()