Eigenproblem for integral operator

On Tue, Jan 18, 2011 at 12:47:22PM -0000, Raphael Kruse wrote:
> New question #141904 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/141904
>
> Hi,
>
> I am looking for a way to implement an integral operator in Python/Dolfin of the following form:
>
> Carleman operator Q : L^2 -> L^2
> Given: an integral kernel q(x,y), where q is symmetric q(x,y) = q(y,x)
>
> For u in L^2 the integral operator Q is given by
> [Qu](x) = \int q(x,y) u(y) dy
>
> In particular, I am interested in determing the eigenvalues and eigenfunctions of Q, i.e.
>
> Qu = \lambda u
>
> Thus, for a finite element approximation I need the matrix
>
> \int \int q(x,y) u(x) v(y) dy dx
>
> where u and v are the trial and test function, respectively.
>
> I am new to Dolfin and sorry if there already is an answer to this question.
>
> Thank you for any help!

I don't think this is possible. As I understand you want 1-dimensional
trial and test functions and integrate that over a 2-dimensional
domain. DOLFIN assumes that everything you integrate has the same
dimension as the domain.

One inefficient way to do this would be to use 2-dimensional trial and
test functions and then sum appropriately over the columns and rows of
the matrix to reduce to 1-dimensional. I haven't thought it through in
detail but it should be possible.

--
Anders

Hi,

thank you for your answer!

Yes, you are right. The final integral is over a 2-dimensional domain, while the trial and testfunctions are 1-dimensional.

But I hoped one can calculate the desired matrix in a perhaps iterative fashion.
Say we have q(x,y) = q(x -y) such that

[Qu](x) = \int q(x - y) u(y) dy

is nothing more but the convolution of q and u.

Is it possible to represent [Qu](x) as function instance in Fenics for a given trial function u such that

L = Qu*v*dx
c = assemble(L)

and c is the column of my desired matrix, which is related to my trial function u?

Thanks a lot!
Raphael

On Wed, Jan 19, 2011 at 10:39:03AM -0000, Raphael Kruse wrote:
> Question #141904 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/141904
>
> Status: Answered => Open
>
> Raphael Kruse is still having a problem:
> Hi,
>
> thank you for your answer!
>
> Yes, you are right. The final integral is over a 2-dimensional domain,
> while the trial and testfunctions are 1-dimensional.
>
> But I hoped one can calculate the desired matrix in a perhaps iterative fashion.
> Say we have q(x,y) = q(x -y) such that
>
> [Qu](x) = \int q(x - y) u(y) dy
>
> is nothing more but the convolution of q and u.
>
> Is it possible to represent [Qu](x) as function instance in Fenics for a
> given trial function u such that
>
> L = Qu*v*dx
> c = assemble(L)
>
> and c is the column of my desired matrix, which is related to my trial
> function u?

I don't see how that would work. Your L (and c) would then be
something that depends on x. We can only assemble numbers (scalars,
vectors, matrices).

Perhaps if you express your q as a sum of basis functions? Then one
could translate the dependency on x into a dependency on an index so
that the above would be assembled into a matrix? Then another trick or
two might solve the problem... I'm not trying to be cryptic, but I
don't have time to think this through in detail atm.

--
Anders

Hi,

I wrote the following code to solve the problem. I would be very grateful if someone with experience could go over it. Very likely it is possible to make the code more efficient.

"""
This scripts builds a FEM-approximation of a Carleman operator Q, which is
given by an integral kernel q(x,y). Then it solves the corresponding
eigenproblem of finding all eigenpairs (lambda,u) such that Qu = lambda*u.

Q is given by [Qu](y) = \int q(x,y) u(x) dx

Author: Raphael Kruse (Univerity Bielefeld), 26 Jan 2011
"""

from dolfin import *
import numpy as np
from numpy import linalg as LA

# parameters for the kernel q
sigma = 1 # standard deviation
gamma = 0.2 # correlation length

dof = 200 # degrees of freedom
tol = 1E-14 # tolerance for eigenvalues

# Defining the integral kernel
q_1d = Expression('sigma*sigma*exp(-pow(x[0]-y,2)/(gamma*gamma))')
q_1d.sigma = sigma
q_1d.gamma = gamma

# Test for PETSc and SLEPc
if not has_la_backend("PETSc"):
    print "DOLFIN has not been configured with PETSc. Exiting."
    exit()

if not has_slepc():
    print "DOLFIN has not been configured with SLEPc. Exiting."
    exit()

# Define mesh, function space
mesh = UnitInterval(dof-1)
coor = mesh.coordinates()

V = FunctionSpace(mesh, "CG", 1)
v = TestFunction(V)
u = TrialFunction(V)

# building the mass matrix M as PETScMatrix

print "Building the mass matrix M...",
M = PETScMatrix()
a = u*v*dx
assemble(a, tensor=M)

print "Done"

# Building the Qu matrix for u running through the testfunctions
print "Building the Q matrix...",
Qu = np.zeros( (dof,dof) ) # holds all values of convolutions of q with
                            # testfunctions
for i in range(dof):
    q_1d.y = coor[i][0] # y runs through the coordinates of the mesh
    L = q_1d*v*dx
    b = assemble(L, mesh=mesh)
    Qu[i,:] = b.array()

Qu = Qu.transpose()

# Computing the values of \int Qu(y) v(y) dy for each u,v

Q = np.zeros( (dof, dof) ) # the matrix Q
Fct_Qu = Function(V)
for k in range(dof):
    quk = np.array(Qu[k][:]) # runs through the values of u
    Fct_Qu.vector()[:] = quk
    L = Fct_Qu*v*dx
    b = assemble(L, mesh=mesh)
    Q[:,k] = b.array()

print "Done"

# transforming Q into PETScMatrix
print "Transforming Q into PETScMatrix format...",
Qh = PETScMatrix(M)

for i in range(dof):
    for k in range(dof):
        Qh[i,k] = Q[i,k]

print "Done"
#print Qh.array()

# Create eigensolver
eigensolver = SLEPcEigenSolver()
eigensolver.parameters["solver"] = "lapack"

# Compute all eigenvalues of Q x = \lambda M x

Part 2 of the code:

# Compute all eigenvalues of Q x = \lambda M x
print "Computing eigenvalues...",
eigensolver.solve(Qh,M)
print "Done"

# Extract largest (first) eigenpair

i = 0
while i < dof:
    r, c, rx, cx = eigensolver.get_eigenpair(i)
    print "%2i. eigenvalue: %g" % (i+1, r)
    i += 1
    if abs(r) < tol:
        i = dof

# Plot eigenfunction corresponding to the largest eigenvalue
r, c, rx, cx = eigensolver.get_eigenpair(2)
print "Eigenvalue corresponding to the plotted eigenfunction is: %g" \
        % r

eigfct = Function(V, rx)
plot(eigfct, title='Eigenfunction')

# Testing the eigenfunction
Qu = np.zeros( dof ) # holds all values of convolutions of q with
                                  # testfunctions
for i in range(dof):
    q_1d.y = coor[i][0] # y runs through the coordinates of the mesh
    L = q_1d*eigfct*dx
    Qu[i] = assemble(L, mesh=mesh)

Qufct = Function(V)
Qufct.vector()[:] = np.array(Qu)

plot(Qufct, title='Qu')

interactive()

This question was expired because it remained in the 'Open' state without activity for the last 15 days.