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.