The error between the exact solution and FEM solution in FEniCS is too large

Question

As indicated in the subject, would some body please indicate the problem in my coding.

What am I missing here?

my code:

""" This solves:

           y''(x) + 4y(x) = 0

           y(0) = -2, y(Pi/4) = 10

    The Exact solution to this problem is:

          y(x) = -2cos(2x) + 10sin(2x))
 """
from dolfin import *
import numpy

a = pi/4.0
mesh = IntervalMesh(60, 0, a)
V = FunctionSpace(mesh, 'CG', 4)

# For left boundary, meaning for y(0) = - 2
ul = Expression("-2- x[0]")
def left_boundary(x, on_boundary):
    tol = 1e-14    
    return on_boundary and abs(x[0])<tol
bc1 = DirichletBC(V, ul, left_boundary)

# For right boundary, meaning for y(pi/4) = 10
ur = Expression("10 - x[0] ")
def right_boundary(x, on_boundary):
    tol = 1e-14    
    return on_boundary and abs(x[0] - pi/4) < tol
bc2 = DirichletBC(V, ur, right_boundary)

# combing bc1 and bc2 in the form of list.
bc = [bc1, bc2]

u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0.0)
a = inner(grad(u), grad(v))*dx + 4.0*u*v*dx
L = f*v*dx

u = Function(V)
solve( a == L, u, bc)

#plot(u, interactive = True)

# Examing the discrete solution
u_nodal_values = u.vector() # the result is a dolfin vector, having all nodal values.
u_array = u_nodal_values.array() # converting the dolfin vector to python array for further processing.

# printing the solution at nodes, meaning u(x_i)
coor = mesh.coordinates()

# comparing the finite element solution with the exact solution.
# finding the maximum error.

# First defining the exact solution.
ue = Expression("-2*cos(2*x[0]) + 10*sin(2.0*x[0])")
u_e = interpolate(ue, V)
u_e_array = u_e.vector().array()
print "Max. error:", numpy.abs(u_e_array - u_array).max()

# now printing the approxmiate solution, the exact solution and error side by side.
#if mesh.num_vertices() == len(u_array):
#   for i in range(mesh.num_vertices()):
#        print 'u(%8g)=%g \t %g \t %8g' %(u(coor[i][0]), u_array[i], u_e_array[i], u_e_array[i]-u_array[i])

Answer 1

Hi, your u_r at $\pi/4$ does not evaluate to 10, so you are prescribing wrong value. The sign of the term in a coming from integration by parts should be negative.

ur = Constant(10)
# ...
a = -inner(grad(u), grad(v))*dx + 4.0*u*v*dx

The simplest way to compare your solution with exact one is via errornorm

print 'L2 norm of error', errornorm(ue, u)

Finally, the last commented out loop might work here but it is a bad practise. You should use dofmaps.

Comments

Thank you for responding.

I am a bit confused here.

You said in my earlier question:

  "for u(0) = 1 use for example u0 = Expression("1-x[0]") "

So this is what I am doing here. Would you please explain it a little. What is the general principle?

Thank you

---

Note the 'for example' in my answer. If you want u(0)=1 then the simplest way is to set
u0=Constant(1) since wherever (any $x$) you evaluate u0 you'll get 1. You can however also set u0 = Expression('1-x[0]'), u0 = Expression('1+x[0]') or u0=Expression('cos(x[0])') since all those expressions evaluate to 1 when $x=0$. That is the princile.

---

So what I have understood is:

You can do ul = Expression( 1- x[0]), for example, since you left boundary is at x = 0. This means here value of 'x[0]' is 0.

and

now since my right boundary is at x = pi/4. ur = Expression(10 - x[0]), would now corresponds to value of u1 at x = pi/4. Further 'x[0]' would now act as a variable whose value is pi/4.

But where am I telling my script to use x[0] as 0 and then pi/4?

Perhaps the statment:

bc1 = DirichletBC( V, ul, boundary) does that for ul.
and
bc2 = DirichletBC( V, ur, boundary) does that for ur.

Correct??

P.S to admins, I am really having hard time typing , since my cursor is hanging a lot. I even changed my browser while writing this.

---

Yes, the third argument to DirichletBC tells it where to evaluate ul and ur. If you

plot(ur, mesh=mesh)

ur is evaluated at many different points to create the plot.

---

Thank you for your time and teaching me some new things!!