Implicit Euler Scheme

Hi. I’m trying to evolve a complex gaussian package under the harmonic oscillator Hamiltonian by means of an semi-implicit Euler scheme. In the next code at final line I received the error message “IndexError: tuple index out of range” but I don’t found that’s happen:

from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as mpl

a = -4.5
b = 4.5
N = 200
T = 6.0            # final time
num_steps = 80     # number of time steps
n_save = 5            # save every n_save steps
Dt = T / num_steps # time step size
x0 = -1.         # initial coordinate
p0 = 0.5        # initial momentum

mesh = generate_mesh(Rectangle(Point(a, -1), Point(b,1)), 64)# mesh has 2 dimensions: real and imaginary parts
V = FunctionSpace(mesh, 'CG', 2)

tol = 1E-14
def boundary(x, on_boundary):
      return on_boundary
bc = DirichletBC(V, Constant(0), boundary)

u_0 = Expression('cos(p0*x[0])*exp(-0.5*pow((x[0]-x0),2))', 'sin(p0*x[0])*exp(-0.5*pow((x[0]-x0),2))',degree=2, p0=p0, x0=x0)
    
u_n = Function(V)
u_n.assign(u_0)

v = TestFunction(V)
u = TrialFunction(V)

p = Expression('x[0]*x[0]',degree=2)
a = (u[0]*v[0]+u[1]*v[1] + Dt/4.*(dot(grad(v[0]),grad(u[1])) - dot(grad(v[1]),grad(u[0]))) + Dt/4. *p*(u[1]*v[0]+u[0]*v[1]))*dx

Several things in this code strikes me as odd:

As you have not encapsulated your code with ``` this might be an issue with formatting. However, to me this looks like you are trying to create a vector expression. However, your function space is a scalar space.
Similarly, this is an issue in your code when you define your variational form as well.
If you want to have a vector function, you need to create a VectorFunctionSpace.
i.e. replace:

with

V = VectorFunctionSpace(mesh, "CG", 2)

This will create a second order vector function space (with an x and y component).
Additionally you need to account for your vector space in the following definitions

bc = DirichletBC(V, Constant((0, 0)), boundary)
u_0 = Expression(("cos(p0*x[0])*exp(-0.5*pow((x[0]-x0), 2))",
                  "sin(p0*x[0])*exp(-0.5*pow((x[0]-x0), 2))"), degree=2, p0=p0, x0=x0)
1 Like

Thanks, this works fine. Other problem was happened. In the last line on my code:

from fenics import *
from mshr import *
import numpy as np
import matplotlib.pyplot as mpl

mesh = generate_mesh(Rectangle(Point(-5, -1), Point(5,1)), 64)
V = VectorFunctionSpace(mesh, "CG", 2)

def boundary(x, on_boundary):
                return on_boundary
bc = DirichletBC(V, Constant((0, 0)), boundary)

u_0 = Expression('cos(p0*x[0])*exp(-0.5*pow((x[0]-x0),2))', 'sin(p0*x[0])*exp(-0.5*pow((x[0]-x0),2))',degree=2, p0=p0, x0=x0, domain=mesh)

u_n = Function(V)
u_n.assign(u_0)

I received the message
*** Error: Unable to interpolate function into function space.
*** Reason: Rank of function (0) does not match rank of function space (1).
*** Where: This error was encountered inside FunctionSpace.cpp.

That is because you Did not copy the last bit of code from my post

There is an extra set of brackets here compared to your code.