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.
-
While Fenics is outputting a solution it seems to be much larger in magnitude to the analytical solution.
-
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()
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
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!