Higher Order Eigenproblem in 1D

General
1

Hi, I am trying to solve a higher order Eigenproblem, i.e. one which contains fourth order derivatives. Specifically, as a simple example I tried to find the Eigenvalues of the following:

\frac{\partial^4 v}{\partial x^4} = \lambda \frac{\partial^2 v}{\partial x^2},\quad v(\pm 1) = v''(\pm 1) = 0

Which should have eigenvalues, \lambda_n = n^2\pi^2 . My attempt to solve this was to introduce a new variable u = v'', and create the weak form from these two equations together. But the answers I am getting are incorrect. Here is a MWE of the code:

from fenics import *
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as la

# Define the mesh
Nmesh = 1001
mesh = IntervalMesh(Nmesh, -1.00, 1.00)

# Define the function space
F1 = FiniteElement('P', "interval", 1)
element = MixedElement([F1, F1])
V = FunctionSpace(mesh, element)

# Set Boundary Conditions for v and u
bcs = [DirichletBC(V.sub(0), Constant((0.0)), 'on_boundary'),\
       DirichletBC(V.sub(1), Constant((0.0)), 'on_boundary')]

p = TrialFunction(V)
f = TestFunction(V)

# Split functions as necessary
v, u = split(p)
f_1, f_2 = split(f)

# Define the bilinear form
a = (-u.dx(0)*f_1.dx(0) +u*f_2+v.dx(0)*f_2.dx(0))*dx
m = (u*f_1)*dx

# Create Matricies
A, M = PETScMatrix(), PETScMatrix()
b = PETScVector()

# Assemble the matricies
assemble(a, tensor=A)
assemble(m, tensor=M)

# Appply boundary conditions
[bc.apply(M) for bc in bcs]                                                                                   

bcinds = []                                                                                    
for bc in bcs:                                                                                 
    bcdict = bc.get_boundary_values()                                                          
    bcinds.extend(bcdict.keys())   

# This just converts PETSc to CSR                                                              
A = sp.csr_matrix(A.mat().getValuesCSR()[::-1])                                                
M = sp.csr_matrix(M.mat().getValuesCSR()[::-1])                                                

# Create shift matrix                                                                          
shift = 1.2345e10*np.ones(len(bcinds))                                                         
S = sp.csr_matrix((shift, (bcinds, bcinds)), shape=A.shape)     

# Compute Eigenvalues
ev, ew = la.eigs(A+S, 6, M, sigma=0.0,which='SM')                                                              

print(ev)
2

Hi, this high-order problem is related to Euler-Bernoulli beam models. Such models require C1 continuity of the discretized field. This is typically done with Hermite elements which FEniCS is lacking. A remedy is to resort to a DG scheme such as described in https://fenicsproject.org/olddocs/dolfin/1.6.0/python/demo/documented/biharmonic/python/documentation.html
which you would have to adapt to the eigenproblem.

Another remedy is to resort to a relaxed formulation which is more or less what you wanted to achieve, except that you cannot enforce the constraint explicitly but have to penalize in the energy. In mechanics, this corresponds to a Timoshenko beam model. The eigenproblem implementation is described here:
https://comet-fenics.readthedocs.io/en/latest/demo/beam_buckling/beam_buckling.html

1 Like