Dear Fenics community,
I would like to solve the following ODE for $u(x)$ with Fenics:
$(K_2(x)u''(x))''-(K_1(x)u'(x))'=0$ on $0\le x\le 1$
with the BC
$u(0)=0$
$u(1)=0.1$
$u''(0)=0$
$u''(1)=0$.
$K_1(x)$ and $K_2(x)$ are given functions. I will drop the argument $x$ now. Writing out the derivatives gives
$2 K_2' u'''+K_2 u'''' + K_2'' u'' - K_1' u' - K_1 u'' = 0 $.
Instead of one fourth order ODE, Fenics can handle two second order ODE. Thus, I decompose into two second order ODE by
$u''=v$,
$2 K_2' v' + K_2 v'' + K_2 '' v - K_1' u' - K_1 v = 0 $.
Now I multiply each ODE with a test function, the first one with $\delta u$ and the second one with $\delta v$,take the integral over the domain, eliminate all second derivatives by integration by parts, where I use that the test functions vanish at the boundaries:
$\int_0^1 -u' \delta u' - v \delta u dx = 0$
$\int_0^1 2 K_2' v' \delta v - K_2 v' \delta v' + K_2 v \delta v -K_1' u' \delta v -K_1 v \delta v dx=0$
From the mixed Poisson example I have seen that these two equations are just added. I don't really understand whether/how this is works. Anyway, after doing so, when I try to implement the weak form, Fenics complains about a dimension mismatch, which comes from the first and fourth term in the second integral equation. Here is a minimum running example:
from dolfin import *
mesh = IntervalMesh(20, 0, 1)
def DirichletBoundary(x, on_boundary):
return on_boundary
u0=Expression('0.1*x[0]')
v0=Expression('0.0')
FS1 = FunctionSpace(mesh, "CG", 4)
FS2 = FunctionSpace(mesh, "CG", 4)
W = FS1 * FS2
bc = [DirichletBC(W.sub(0), u0, DirichletBoundary),DirichletBC(W.sub(1), v0, DirichletBoundary)]
(u,v) = TrialFunctions(W)
(du,dv) = TestFunctions(W)
K1=Expression('100.0*(x[0]+1.0)')
K2=Expression('500.0*(x[0]+1.0)')
dK1=Expression('100.0')
dK2=Expression('500.0')
ddK1=Expression('0.0')
ddK2=Expression('0.0')
a= (
-inner(grad(u),grad(du))
-inner(v,du)
#~ -2.0*dK2*v*grad(v)[0]
-K2*inner(grad(v),grad(dv))
+ddK2*inner(v,dv)
#~ -dK1*grad(u)[0]*dv
-K1*inner(v,dv)
)*dx
f = Constant(0.0)
L = f*dv*dx
w = Function(W)
solve( a == L, w, bc)
(u, v) = w.split()
plot(u, interactive = True)
In its present form, the code works, but I need to include the missing lines in the definition of a.
The inner Function does not work due to a dimension mismatch that occurs actually only if the problem is higher dimensional. The lines
#~ -dK1*grad(u)[0]*dv
give problems at runtime. I also tried to use derivative and norm to make grad(u) scalar, without success.
Can anyone suggest how to solve this?
