Large Error in Numerical Solution and Confusion about Test Function Dependancy

variational formulation
1

Hi there. I picked Fenics up a couple of days ago and have already managed to get some help from the community so for that I’d like to say thanks a lot!

I am simulating a pair of first order coupled, linear ODEs on the interval [0,\pi]

My ODE’s are:
u_1'=2u_2 and u_2'=-2u_1 subject to u_1(0)=u_1(\pi)=0 and u_2(0)=u_2(\pi)=1. The analytic solution to these equations are u_1=\sin(2x), u_1=\cos(2x).

I have two problems with my implementation.

  1. While Fenics is outputting a solution it seems to be much larger in magnitude to the analytical solution.

  2. As far as I can see there are two ways to form the variational problem with a test function v=(v_1,v_2). That is multiply u_1'=2u_2 by v_1 and u_2'=-2u_1 by v_2 or vice versa, sum and then integrate, both which, as far as I can understand, seem to pose the same variational problem theoretically.

\int v_1(u_1'-2u_2)+v_2(u_2'+2u_1)dx=0, \forall v=(v_1,v_2)
or
\int v_2(u_1'-2u_2)+v_1(u_2'+2u_1)dx=0, \forall v=(v_1,v_2)

However give different solutions when implemented into the code (see below).

My attempt to implement the numerical solution is given below:

from fenics import *
import matplotlib.pyplot as plt
import numpy as np
# Create mesh and define function space
mesh = IntervalMesh(50,0,pi)
#we are solving for two unknowns u1, u2 each ui must be in its own funtion space

P1 = FiniteElement('P',interval, 1)
element = MixedElement([P1,P1])

V = FunctionSpace(mesh,element)

#need to create test and trial functions
v_1,v_2 = TestFunctions(V)

u = TrialFunction(V)
u_1,u_2=split(u)

#defines boundaries
top='near(x[0],pi)'
bottom='near(x[0],0)'

#sets bc at the boundaries
bc_top=DirichletBC(V,[Constant(0),Constant(1)],top)
bc_bottom=DirichletBC(V,[Constant(0),Constant(1)],bottom)
bc_list=[bc_top,bc_bottom]

#now define the variational problem
F=((u_1.dx(0)-2*u_2)*v_2+(u_2.dx(0)+2*u_1)*v_1)*dx

a1=lhs(F)
L1=Constant(0)*v_2*dx 
#need to define manually since Fenics rhs will not assign zero value form

A1=assemble(a1)
b1=assemble(L1)
[bc.apply(b1) for bc in bc_list]

#now bilinear form is defined turn u into a function
u = Function(V)
u_1,u_2=split(u)

solve(A1,u.vector(),b1)

Problem 1:

Here is the plot and the code plotting them.

u_1,u_2=split(u)
x=np.linspace(0,np.pi-0.001)#to avoid plotting out of domain
u1=np.array([u(point)[0] for point in x])
u2=np.array([u(point)[1] for point in x])
plt.grid(True)
plt.plot(x,u1)
plt.plot(x,u2)
plt.show()

gra2

As you can see the numerical solutions have a much higher amplitude than the analytic solutions and what’s more u_2 doesn’t obey the boundary condition u_2(0)=u_2(\pi)=1. I am at a loss as to what’s causing this.

Problem 2:

As mentioned above I am confused which variational formula to input into Fenics (given in code):

A: F=((u_1.dx(0)-2*u_2)*v_1+(u_2.dx(0)+2*u_1)*v_2)*dx
or
B: F=((u_1.dx(0)-2*u_2)*v_2+(u_2.dx(0)+2*u_1)*v_1)*dx

Where v=(v1,v2) is defined in the code via
v_1,v_2 = TestFunctions(V)
The plots given by both of these variational formulations are
gra1gra2

          A                      B

I can see that plot B is the one which agrees most with the analytical solution but I cannot understand why formulation B seems to be the correct variational formulation of the problem to implement within Fenics.

Thanks a lot for anyone who is willing to give some input into my problem!

2

Hi, your solution is wrong because you are not apply boundary conditions to the bilinear form. Consider either:

A1=assemble(a1)
b1=assemble(L1)
for bc in bc_list:
    bc.apply(A1)
    bc.apply(b1)

or even better

A1, b1 = assemble_system(a1, L1, bc_list)

With the proper solution, the choice of the TestFunctions order does not matter