0

I've looked all over stack exchange and can see the solution to this so here's hopping somebody knows. I'm Trying to solve an time-dependent parabolic equation using mixed finite element method and backwards Euler. So I have two coupled equations:

${\partial u \over \partial x} + \alpha \omega = \gamma{\partial \omega \over \partial t} $

${\partial \omega \over \partial x} + \beta u = 0$

Where $u$ and $\omega$ are the parameters to be solved for and $\alpha$, $\beta$ and $\gamma$ are constant (not depending on either $u$ or $\omega$). We can reduce the system down to a matrix equation of the form:

$M_1 \psi^t = M_2 \psi^{t-1}$

where $\psi = (u\:\:\omega)^\text{T}$ and solve like normal. However this is not stably for every $\Delta t$ (time stepping) - only a small range of $\Delta t$ give the correct result i.e if $t_{min} <\Delta t < t_{max} $ . $t_{max}/t_{min}$ obviously this is varying with $\Delta x$, $\alpha$, $\beta$ and $\gamma$.

So my question is:

1) I've double checked my algorithm and believe it to be correct so why is backwards Euler not always stable for some values of $\Delta t$? How does backwards Euler deal with mixed finite elements? Is it still stable?! What am I missing?

In mixed finite elements there are two conditions that must be met: 1) the inf-sup condition 2) the invertibility on the kernel condition

I believe both are still intact - however these are derived for a stationary state (no time dependence) so maybe they are different for time dependence?

I'm really struggling on this so any help is really appreciated. I've been on this for nearly 4 days and I must be missing something. As far as I know the Courant-friedrichs-levy condition doesn't apply to backwards Euler. Also it should be noted that setting $\gamma = 0$ (non-time dependent) gives the correct solutions - so the code is working-ish. Many thanks guys!

user2350366
  • 513

1 Answers1

0

Why are you writing this as two equations?

You could just substitute the second one into the first, and get $$ \frac{\partial{\omega}}{\partial t} = \frac{\alpha}{\gamma}\omega -\frac{\beta}{\gamma} \frac{\partial^2 \omega}{\partial^2 x} $$

This is a standard diffusion equation with source term. There is a constrain for the equation to be stable: $$ -\frac{\beta}{\gamma} > 0 $$ otherwise, this is anti-diffusion, which is a ill-posed problem. Such a ill-posed problem would always blow up, even if you're using an unconditionally stable method like backward-Euler. It blows up because the equation blows up.

Further, watch your boundary conditions. If you want to have a stationary solution, then you cannot use insulated boundaries, i.e. $\frac{\partial \omega}{\partial x} = 0$ on both sides.

rcomblen
  • 166