Solving 1 space 1 time equation with Neumann bc vs Dirichlet bc

1

Hi,

I have been trying to solve the 1-space 1-time problem u_t=u_{xx} on [0,\infty)\times[0,2\pi] subject to the Neumann conditions u'(t,0)=u'(t,2\pi)=e^{-t} and u(0,x)=\sin(x). This has analytical solution u(t,x)=e^{-t}\sin(x) which looks like a sine curve being progressively damped.

My code runs and plots the following curves for t\in\{0,0.1,0.2,...,1\}

download

Which isn’t what the analytical solution looks like. Below is the code used to solve the problem.

from fenics import *
from dolfin import*
import matplotlib.pyplot as plt
import numpy as np

Intervals=10
Time=1
dt=Time/Intervals
t=0

#defines the mesh - the interval
upper=2*np.pi
mesh=IntervalMesh(100,0,upper)
V=FunctionSpace(mesh,'P',1)

#our boundary is 0-dimensional, hence neumann bc have to be manually implemented 
#since boundary normal direction is not well defined!


boundaries = MeshFunction('size_t', mesh,1)


#defines boundary classes
class BoundaryRight(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (near(x[0],upper))

class BoundaryLeft(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and (near(x[0],0))

#creates instances to use in the program
boundaryLeft = BoundaryLeft()
boundaryRight = BoundaryRight()

#marks the boundaries 
boundaryLeft.mark(boundaries,1)
boundaryRight.mark(boundaries,2)
ds = ds(subdomain_data=boundaries)



#give initial t=0 cond and neumann bc.
u0=Expression('sin(x[0])',degree=1)
g=Expression('exp(-t)',degree=1,t=t)
#defines initial u_n as u0 as a function on the mesh
u_n=interpolate(u0,V) 


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

#define the variational problem

a=(u*v+dt*u.dx(0)*v.dx(0))*dx
#minus at ds(1) indicates normal points to the left
L=dt*(+g*v*ds(2) - g*v*ds(1))+u_n*v*dx

for n in range(Intervals):
    t+=dt
    g.t=t
    #solves
    solve(a==L, u)
    #reassigns u_n for next loop iteration
    u_n.assign(u)

Now after some thought it does seem to be the case that my Neumann boundary conditions are being implemented, the solution is becoming flatter at the ends of the interval as t \to \infty. However the solution is off. I tested the same code but instead with the Dirichlet boundary conditions u(t,0)=u(t,2\pi)=0 and hey presto:
download (1)

Which does agree with the analytical solution. So I am wondering:

  1. Am I implementing the Neumann boundary conditions wrong?

or

  1. Do I need to specify more bc’s in the Neumann problem for the problem to be well posed?

In both the Neumann and Dirichlet problems I have specified 3 boundary conditions, the initial curve at t=0 and the bc’s at the edges of the interval. However when these last two bc’s have been Dirichlet I have gotten a unique solution and when they have been Neumann my solution doesn’t make a lot of sense. Do Dirichlet bc’s carry more information than Neumann bc’s? Or maybe even Fenics needs more Neumann bc’s than it would Dirichlet?

Thank you for any feed back you are willing to offer!

2

The mistake is in your definition of the mesh function.
As you are constraining points, you need to use the following mesh function:

boundaries = MeshFunction('size_t', mesh,0)

as

boundaries = MeshFunction('size_t', mesh,1)

is the function marking individual cells.
Using the correct marker solves your issue

1 Like