### Solving u_tt - u_xxxx = f(t,x) by modifying biharmonic demo.

83

views

0

I'm trying to solve the following PDE:

$u_{tt}-u_{xxxx}=f\left(t,x\right)$

with boundary conditions in space:

$u\left(t,0\right)=u\left(t,1\right)=0$

$u_{xx}\left(t,0\right)=u_{xx}\left(t,1\right)=0$

and Neumann boundary conditions in time:

$u_t\left(0,x\right)=g_0\left(x\right)$

I'm trying to modify the biharmonic demo, but don't know how to properly treat the time variable.

I'm not an expert in FEM, my question may lack of basic knowledges, sorry. Any advice and suggestion are welcome. Thanks.

Here's what I got so far:

$u_{tt}-u_{xxxx}=f\left(t,x\right)$

`u`_{tt}−`u`_{xxxx}=`ƒ`(`t`,`x`)with boundary conditions in space:

$u\left(t,0\right)=u\left(t,1\right)=0$

`u`(`t`,0)=`u`(`t`,1)=0$u_{xx}\left(t,0\right)=u_{xx}\left(t,1\right)=0$

`u`_{xx}(`t`,0)=`u`_{xx}(`t`,1)=0and Neumann boundary conditions in time:

$u_t\left(0,x\right)=g_0\left(x\right)$

`u`_{t}(0,`x`)=`g`_{0}(`x`) , $u_t\left(1,x\right)=g_1\left(x\right)$`u`_{t}(1,`x`)=`g`_{1}(`x`)I'm trying to modify the biharmonic demo, but don't know how to properly treat the time variable.

I'm not an expert in FEM, my question may lack of basic knowledges, sorry. Any advice and suggestion are welcome. Thanks.

Here's what I got so far:

```
from fenics import *
import numpy as np
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True
# Make mesh ghosted for evaluation of DG terms
parameters["ghost_mode"] = "shared_facet"
# Create mesh and define function space
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "CG", 2)
# Define Dirichlet boundary
def boundary(x):
return near(x[1], 0.0) or near(x[1], 1.0)
# Define boundary condition
bc = DirichletBC(V, 0.0, DirichletBoundary())
# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
# Define normal component, mesh size and right-hand side
h = CellDiameter(mesh)
h_avg = (h('+') + h('-'))/2.0
n = FacetNormal(mesh)
f = Expression('-(pi*pi+pi*pi*pi*pi)*sin(pi*x[0])*sin(pi*x[1])', degree=2)
g = Expression('pi*sin(pi*x[0])*cos(pi*x[1])', degree=2)
# Penalty parameter
alpha = Constant(8.0)
# Define bilinear form
a = ( inner(u.dx(0),v.dx(0)) + inner(u.dx(1).dx(1), v.dx(1).dx(1)) )*dx \
- inner(avg(div(grad(u))), jump(grad(v), n))*dS(1) \
- inner(jump(grad(u), n), avg(div(grad(v))))*dS(1) \
+ alpha/h_avg*inner(jump(grad(u),n), jump(grad(v),n))*dS(1)
# Define linear form
L = -f*v*dx + g*v*ds(0)
# Solve variational problem
u = Function(V)
solve(a == L, u, bc)
```

Community: FEniCS Project

### 1 Answer

0

You can always discretize the second-order time derivative with one of the finite difference methods:

- backward difference:

\[ \frac{\partial^2 u^n}{\partial t^2} \approx \frac{\partial}{\partial t} \left( \frac{u^n - u^{n-1}}{\Delta t} \right) = \frac{1}{\Delta t} \left( \frac{\partial u^n}{\partial t} - \frac{\partial u^{n-1}}{\partial t} \right) \approx \frac{1}{\Delta t} \left( \frac{u^n - u^{n-1}}{\Delta t} - \frac{u^{n-1} - u^{n-2}}{\Delta t} \right) = \frac{u^n - 2 u^{n-1} + u^{n-2}}{(\Delta t)^2} \] - forward difference:

\[ \frac{\partial^2 u^n}{\partial t^2} \approx \frac{\partial}{\partial t} \left( \frac{u^{n+1} - u^n}{\Delta t} \right) = \frac{1}{\Delta t} \left( \frac{\partial u^{n+1}}{\partial t} - \frac{\partial u^n}{\partial t} \right) \approx \frac{1}{\Delta t} \left( \frac{u^{n+2} - u^{n+1}}{\Delta t} - \frac{u^{n+1} - u^n}{\Delta t} \right) = \frac{u^{n+2} - 2 u^{n+1} + u^n}{(\Delta t)^2} \] - central difference:

\[ \frac{\partial^2 u^n}{\partial t^2} \approx \frac{\partial}{\partial t} \left( \frac{u^{n+1/2} - u^{n-1/2}}{\Delta t} \right) = \frac{1}{\Delta t} \left( \frac{\partial u^{n+1/2}}{\partial t} - \frac{\partial u^{n-1/2}}{\partial t} \right) \approx \frac{1}{\Delta t} \left( \frac{u^{n+1} - u^n}{\Delta t} - \frac{u^n - u^{n-1}}{\Delta t} \right) = \frac{u^{n+1} - 2 u^n + u^{n-1}}{(\Delta t)^2} \]

Please login to add an answer/comment or follow this question.

Thanks for your reply. Considering the Neumann boundary conditions on two sides, it seems not easy to apply FDM scheme in this case.

I'm currently trying to decouple the PDE into two heat equations.