FEniCS incorrectly evaluates projection

[Glen D. Granzow]

The following lines should create a function g which is the projection of
f(x) = x^2 on a one-dimensional function space with Lagrange elements of
order one. I do not think that FEniCS does this correctly.

##### Beginning of python code #####

from dolfin import *

n = int(raw_input('Please enter a value for n: '))

mesh = UnitIntervalMesh(n)
f = Expression('x[0]*x[0]')
V = FunctionSpace(mesh, 'Lagrange', 1)

u = TrialFunction(V)
v = TestFunction(V)
g = Function(V)
solve(u*v*dx == f*v*dx, g)
print 'g.vector().array() =', g.vector().array()

# The same result can be obtained from

g = project(f,V)
print 'g.vector().array() =', g.vector().array()

##### End of python code #####

If n = 1, for example, the output produced is

g.vector().array() = [ 0. 1.]

indicating that the function g is

g(x) = 0*v_0(x) + 1*v_1(x) = 0*(1-x) + 1*x = x

which linearly interpolates between f(x) evaluated at the nodes, f(0)=0 and
f(1)=1. I think that the projection of f(x) onto V should be

g(x) = (-1/6)*v_0(x) + (5/6)*v_1(x) = (-1/6)*(1-x) + (5/6)*x = x - 1/6

HERE IS MY REASONING: The goal is to find an approximate solution to

g(x) = f(x)

which is written in the weak form

integral(g(x)v(x)) = integral(f(x)v(x)) for all v(x)

When discretized, g is approximated by

G(x) = g_0*v_0(x) + g_1*v_1(x)

The coefficients g_0 and g_0 are determined by solving the system of equations

integral(G(x)*v_0(x)) = integral(f(x)*v_0(x))
integral(G(x)*v_1(x)) = integral(f(x)*v_1(x))

When I do these integrals I get the system of equations

(1/3)*g_0 + (1/6)*g_1 = 1/12
(1/6)*g_0 + (1/3)*g_1 = 1/4

whose solution is

g_0 = -1/6
g_1 = 5/6

CONCLUSION: The projection of f(x) = x^2 on V is G(x) = (-1/6)*v_0 + (5/6)*v_1

For n > 1, FEniCS also returns a function g(x) which linearly interpolates
between f(x) evaluated at the nodes. By my reasoning this is not correct since

G(x) > f(x) for all x in (0,1)\{x_0, x_1, ..., x_n}

(where the x_i's and the nodes) so

integral(G(x)*v_i(x)) > integral(f(x)*v_i(x)) for i = 0,1,2,...,n

###########################################################################

If you disagree with my reasoning or analysis please point out the error of
my ways. Thank you.

[Martin Sandve Alnæs (martinal)]

FEniCS has moved to bitbucket, please see http://fenicsproject.org/support/ for how to ask questions now.

You need to set the degree of f:

    f = Expression('x[0]*x[0]', degree=2)

[Glen D. Granzow (ggranzow)]

Thank you Martin for your prompt attention to my bug report. Using "degree=2" does indeed produce the projection that I expected. Is there somewhere in the FEniCS documentation where the "degree" keyword argument is described? I had no idea that it existed and it now seems important.

[Martin Sandve Alnæs (martinal)]

Should be in the Expression doc string.

Martin
1. okt. 2013 20:35 skrev "Glen D. Granzow" <email address hidden> følgende:

> Thank you Martin for your prompt attention to my bug report. Using
> "degree=2" does indeed produce the projection that I expected. Is there
> somewhere in the FEniCS documentation where the "degree" keyword
> argument is described? I had no idea that it existed and it now seems
> important.
>
> --
> You received this bug notification because you are subscribed to FEniCS
> Project.
> https://bugs.launchpad.net/bugs/1232270
>
> Title:
> FEniCS incorrectly evaluates projection
>
> Status in DOLFIN:
> Invalid
>
> Bug description:
> The following lines should create a function g which is the projection of
> f(x) = x^2 on a one-dimensional function space with Lagrange elements of
> order one. I do not think that FEniCS does this correctly.
>
> ##### Beginning of python code #####
>
> from dolfin import *
>
> n = int(raw_input('Please enter a value for n: '))
>
> mesh = UnitIntervalMesh(n)
> f = Expression('x[0]*x[0]')
> V = FunctionSpace(mesh, 'Lagrange', 1)
>
> u = TrialFunction(V)
> v = TestFunction(V)
> g = Function(V)
> solve(u*v*dx == f*v*dx, g)
> print 'g.vector().array() =', g.vector().array()
>
> # The same result can be obtained from
>
> g = project(f,V)
> print 'g.vector().array() =', g.vector().array()
>
> ##### End of python code #####
>
> If n = 1, for example, the output produced is
>
> g.vector().array() = [ 0. 1.]
>
> indicating that the function g is
>
> g(x) = 0*v_0(x) + 1*v_1(x) = 0*(1-x) + 1*x = x
>
> which linearly interpolates between f(x) evaluated at the nodes, f(0)=0
> and
> f(1)=1. I think that the projection of f(x) onto V should be
>
> g(x) = (-1/6)*v_0(x) + (5/6)*v_1(x) = (-1/6)*(1-x) + (5/6)*x = x - 1/6
>
> HERE IS MY REASONING: The goal is to find an approximate solution to
>
> g(x) = f(x)
>
> which is written in the weak form
>
> integral(g(x)v(x)) = integral(f(x)v(x)) for all v(x)
>
> When discretized, g is approximated by
>
> G(x) = g_0*v_0(x) + g_1*v_1(x)
>
> The coefficients g_0 and g_0 are determined by solving the system of
> equations
>
> integral(G(x)*v_0(x)) = integral(f(x)*v_0(x))
> integral(G(x)*v_1(x)) = integral(f(x)*v_1(x))
>
> When I do these integrals I get the system of equations
>
> (1/3)*g_0 + (1/6)*g_1 = 1/12
> (1/6)*g_0 + (1/3)*g_1 = 1/4
>
> whose solution is
>
> g_0 = -1/6
> g_1 = 5/6
>
> CONCLUSION: The projection of f(x) = x^2 on V is G(x) = (-1/6)*v_0 +
> (5/6)*v_1
>
> For n > 1, FEniCS also returns a function g(x) which linearly
> interpolates
> between f(x) evaluated at the nodes. By my reasoning this is not
> correct since
>
> G(x) > f(x) for all x in (0,1)\{x_0, x_1, ..., x_n}
>
> (where the x_i's and the nodes) so
>
> integral(G(x)*v_i(x)) > integral(f(x)*v_i(x)) for i = 0,1,2,...,n
>
>
> ###########################################################################
>
> If you disagree with my reasoning or analysis please point out the error
> of
> my ways. Thank you.
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/dolfin/+bu...

[Glen D. Granzow (ggranzow)]

Yes, I see that the "degree" argument is mentioned in the Expression doc string. But there is no discussion of its usage, default, or significance. I was really taken by surprise by the necessity of specifying this since after reading the tutorial I had the impression that FEniCS was going to integrate expressions exactly (that is, symbolically). Perhaps this means than an expression like "cos(x[0])" will never be integrated exactly since an infinite order polynomial is required. Thank you again Martin.