Cannot find bug in very easy fourth order ODE

1

Hi!
I want to solve a very easy fourth order problem with Robin BC on the 1D interval (0,1) in FEniCS. The weak form says:
<u’‘,v’‘> + 1/b_1 * (u(0)v(0)+u(1)v(1)) + 1/b_2*(u’(0)v’(0)+u’(1)v’(1)) = <f,v>

For the rhs f specified in the code below, I know the true solution true_sol. Both Expressions should be right. Nevertheless, when I try to plot the numerical solution in comparison to the true solution, the numerical solution is completely false. All the code is attached below.
As I am a complete FEniCS beginner, I have no clue what I am doing wrong. Probably some very dumb rookie mistake… Thanks very much in advance for any help.

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

n = 500
mesh = IntervalMesh(n, 0, 1)
V = FunctionSpace(mesh, 'CG', 2) 

b1 = 1
b2 = 1
u = TrialFunction(V)
v = TestFunction(V)
f = Expression('pow((2*pi),4)*sin(2*pi*x[0])', degree = 2)


a = dot(div(grad(u)), div(grad(v)))*dx + pow(b1,-1)*(u*v*ds) + pow(b2,-1)*(dot(grad(u),grad(v))*ds)
L = f*v*dx

sol = Function(V)

solve(a==L, sol)

true_sol = Expression('sin(2*pi*x[0]) - 2*pi*(2*x[0]-1)*(4*b1*(pow(pi,2)*(6*b2-2*pow(x[0],2)+2*x[0]+1)+3)+x[0]*(x[0]-1))*pow(24*b1 + 6*b2 + 1, -1)', domain=mesh, degree = 4, b1 = b1, b2 = b2)

#Plots
dom = np.arange(0,1,0.01)
values_sol = []
values_true = []
for i in dom:
    values_sol.append(sol(i))
    values_true.append(true_sol(i))
plt.plot(dom, values_sol, c='r', label = 'Calculated')
plt.plot(dom, values_true,'b', label = 'True')
plt.legend()
plt.savefig('Figure')
2

See for instance: Biharmonic equation — DOLFIN documentation

2 Likes
3

Thanks for the quick answer. Sadly, that file was exactly my starting point: I tried to adjust this to 1D and to my bilinear form/boundary conditions and arrived at the non-working solution above.

4

But you have removed all coupling conditions from your form. (i.e. all dS integrals, those over internal facets).

2 Likes
5

Thanks, don’t know how I could forget that…