Time dependend solving of equations

1

Hi everyone!

I am still quite new to FEniCS and for my Masters Thesis I am wokring with the biharmonic equation. I figured out, that due to the restriction of available finite elements i had to reduce the equations order and that resultated in a system of equations. So far so good, i could already calculate Eigenfunction and somewhat verify them. Now i want to create a standing wave. To solve the system in regards to time I am using Störmers Method, as used for the wave equations in the book “Finite Difference Computing with Python” by Hans Petter Langtangen. Which replaces time derivatives of second order with the following approximation:
(u^n+1 - 2u^n + u-1)/(delta t^2)
and for the first step:
(u^1-u^-1)/(2*delta t)
The problem that now arises is that the solution somehow swings out to zero.
Can somebody can give me a hint or an Ansatz to solve my problem?
Thanks in advance and best regards.
The code is the following:

from __future__ import print_function
from fenics import *

mesh = UnitSquareMesh(20,20)
V = FunctionSpace(mesh,'CG', 2)

# Define Dirichlet boundary
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary

# Define boundary condition
u_0 = Constant(0)
bc = DirichletBC(V, u_0, DirichletBoundary())


# Define trial and test function
u = TrialFunction(V)
v = TrialFunction(V)
phi_1 = TestFunction(V)
phi_2 = TestFunction(V)

W = as_matrix([[v],[u]])
Phi = as_matrix([[phi_1],[phi_2]])

a =  (inner(grad(W), grad(Phi)))*dx 
A = PETScMatrix()
assemble(a, tensor=A)

eigensolver = SLEPcEigenSolver(A)

# Compute all eigenvalues of A x = \lambda x
print("Computing eigenvalues. This can take a minute.")
eigensolver.solve()


# extract eigenvalues and assemble them in a plot
r, c, rx, cx = eigensolver.get_eigenpair(4)
        
print("Largest eigenvalue: ", r)
   
# Initialize function and assign eigenvector
u_lambda = Function(V)
u_lambda.vector()[:] = rx
plot(u_lambda) 
# %% Computation of time dependent problem
import math 

T = 1.5         # final time
num_steps = 300     # number of time steps
dt = T / num_steps # time step size

u_m = Function(V)
u_n = Function(V) 


u = TrialFunction(V)
phi = TestFunction(V)
# Create VTK file for saving solution
vtkfile = File('solution/ResonatingPlate.pvd')


#define F and assemble unknown and known variables
F = (u*phi - 2*u_n * phi + u_m * phi)*dx + dt*dt *a
a_evo, L = lhs(F), rhs(F)

u = Function(V)

# Time-stepping
u_m.assign(u_lambda)
u_n.assign(u_lambda - 0.5*dt*dt*u_lambda)
t = 0
for n in range(num_steps):   
    
    # Compute solution
    solve(a_evo == L, u, bc)
    # Save to file and plot solution
    vtkfile << (u_m, t)
    u_m.assign(u_n)
    # Update previous solution
    u_n.assign(u)
    # Update current time
    t += dt