Crank Nicolson inconsistent with backward Euler, modified heat equation

From what I understood of the FEniCS tutorial, when the heat equation is presented, it is solved with a backward Euler discretization, because even though it is slow, it should converge towards the solution without any problem. Crank Nicolson discretization is also mentioned and implemented later, although its stability is not ensured.

However when I compare the solution to a heat equation (not in the tutorial), they do not converge towards the same solution according to the time discretization I choose. Physically there is a single solution, and as far as I can tell the backward Euler method is the correct one. But I am not sure why the Crank Nicolson method “fails” (it is not obvious at all since the code converges).

The equation to solve is the 1-dimensional version of \nabla \cdot (\kappa \nabla T) -\mu\vec J \cdot \vec \nabla T+\rho J^2 =C_p\frac{\partial T}{\partial t}, where \kappa=\kappa_0+\kappa_1 T and \rho(T) is also linear with respect to temperature, \mu and C_p are constant and \vec J is sinusoidal in time.

Then at each time step I solve for a voltage equation (Ohm’s law in the material), and I stop to compute when the voltage oscillations do not change much anymore. Note that the temperature oscillations caused by the time oscillating value of \vec J have an impact on the voltage oscillations.

The voltages I get are quite different according to if I employ the backward Euler method or the Crank Nicolson one.

My weak forms are:
F_Crank_Nicolson = -C_p*v*(Temp-T_n)*dx + 0.5*( -(kappa_0+kappa_1*Temp)*dt*dot(grad(Temp), grad(v))*dx - mu*dt*J*Temp.dx(0)*v*dx + J**2*(rho_func)*dt*v*dx -(kappa_0+kappa_1*T_n)*dt*dot(grad(T_n), grad(v))*dx - mu*dt*J_n*T_n.dx(0)*v*dx + J_n**2*(rho_func_previous)*dt*v*dx)

and

F_backward_Euler = -(kappa_0+kappa_1*Temp)*dt*dot(grad(Temp), grad(v))*dx - mu*dt*J*Temp.dx(0)*v*dx + J**2*(rho_func)*dt*v*dx -C_p*v*(Temp-T_n)*dx

Between the 2 codes everything else is kept the same, i.e. the mesh, the time step, etc. I have also varied the time step towards smaller and smaller values to see that both methods would converge. They do converge, but to different values. I am wondering this behavior is to be expected, or whether I made a mistake in the weak form of the Crank Nicolson discretization.

If you want to learn more about the stability of temporal discretization schemes, I would suggest reading chapter 1 and 2 of http://hplgit.github.io/INF5620/doc/pub/H14/diffu/pdf/main_diffu-4print.pdf
and chapter 19.10 of http://hplgit.github.io/INF5620/doc/pub/H14/fem/pdf/main_fem-4print.pdf

I would also suggest starting with a simplified version of the heat equation, that can be verified against an analytical solution, to ensure that your discretization is correct, as done in the fenics tutorial: Solving PDEs in Python - <br> The FEniCS Tutorial Volume I

1 Like

Thank you dokken for the references. I have quickly taken a look and it matches what I had read before, i.e. that Crank-Nicolson could be unstable, under certain conditions. But it isn’t clear to me what would happen if the code is unstable. For example, in my case when I pick smaller and smaller time steps, the values I get seem to converge towards some values. It does look like it is stable. What is strange is that these values are quite different from the ones I get with the backward Euler method.

I don’t really see the point of solving the simpler heat equation and compare it with the analytical formula. Sure, it will work, but then what? I am considering a not so complicated equation either, only 2 “simple” terms are added compared to the one of the tutorial.

See for instance chapter 2.8 of: http://hplgit.github.io/INF5620/doc/pub/H14/diffu/pdf/main_diffu-4print.pdf
which I posted in the previous post

1 Like