I'm looking at the parabolic equation
$$
\frac{\partial u}{\partial t} – \frac{\partial}{\partial x}\left(\alpha(u)\frac{\partial u}{\partial x}\right) = f(x,t)
$$
I know the Backwards Euler method is A-stable, meaning that no matter what $\Delta t$ I choose, I should get a stable solution, correct? I'm trying to implement a code to solve the above problem (using a finite element method in space– DG method to be exact, and BE in time) and my issue is that the solution is most definitely not stable, even for small time steps like $\Delta t = 0.01$. I have to make it as small as $\Delta t = 0.000001$ in some cases just to get a decent solution, which is garbage. Anyway, I'm trying to test my code on simpler cases:
1) Linear ($\alpha(u) = 1$), non-time dependent problem with $\Delta t = 100$ and the final time being $t_{f} = 5\Delta t$. No problems here. Solution is stable and converges for mesh refinement study.
2) Non-linear ($\alpha(u) = u$), non-time dependent problem with the same parameters as above. I'm using Newton's method to solve the resulting system. Again, no problems. Solution converges as well.
3) Linear ($\alpha(u) = 1$), time-dependent problem with the same parameters as above. Again, no problems.
4) Non-linear ($\alpha(u) = u$), time-dependent problem with the same parameters as above. VERY BIG PROBLEMS HERE. I need to lower the time step to about 0.01 just to see a decent result.
I've been staring at my algorithm for hours wondering if I did something wrong but now I'm trying to think about the method itself.
Best Answer
Short Answer: Backward Euler method is not A-stable for nonlinear diffusion equation where the diffusion constant is "strongly" dependent on the solution $u$ itself.
Long Answer (but rather heuristic): Backward Euler method (BEM) is stable for $$ \frac{\partial u}{\partial t} = Au + f $$ when the differential operator $A$ has all eigenvalue being negative (or all eigenvalue's real parts are negative). An intuitive explanation is that, the error introduced by the scheme itself in current time step will not propagate to the next time step in a magnified fashion. But rather, the error at $n$-th step will shrink as we marching to the next time step, and the sum of the error in each time adding up will be order $O(h^{k})$.
Your equation is: $$ \frac{\partial u}{\partial t} = \nabla\cdot(\alpha(u)\nabla u) + f $$ written in dimensionless form.
If $\alpha$ is a constant, or is something not relying on $u$, and has certain smoothness. Assuming you have the Dirichlet boundary condition, then the following eigenvalue problem $$ A u = \nabla\cdot(\alpha \nabla u) = \lambda u$$ has all eigenvalues, $\lambda$, are negative. Using BEM will converge nicely.
If $\alpha =\alpha(u)$, the operator you have here is nonlinear $A(u)$. Rewriting the equation as: $$ \frac{\partial u}{\partial t} = [A(u)]u + f = Au + [A(u)-A]u + f $$ where $Au = \nabla\cdot(\hat{\alpha}(x) \nabla u) $. If $\|A(u) - A \|$ is big, then BEM won't be stable, while iterative method works.
Lastly, I suggest you try semi-implicit Runge-Kutta scheme for time dimension. At time step $t_n$, construct an approximation to the operator $A(u({t_n},x))$ using $A(u({t_{n-1}},x))$, $A(u({t_{n-2}},x))$, and so on.