Debugging Newton-method used in a CG-approach

Question

I am currently proof-checking my program, which is intended to use Newton's method for solving nonlinear equations, using a continuous galerkin approach. Thus, as first step I checked it using a time-independent equation ($H^1$-convergence rate of finite element method for Poisson equation, depending on element order), with positive results. As next step I intended to check the time-dependent heat equation $$\partial_tu=\nabla^2u\text{, }u_0=\sin(\pi x)\sin(\pi y)$$ on the domain $$\Omega=[0,1]\times[0,1]$$ and the boundary conditions $$u=0\text{ on }\partial\Omega$$ As time discretization I use $$\frac{u_{new}-u_{old}}{dt}=\nabla^2u_{new}$$ This should give me $$u(x,y,t)=\exp(-2\pi^2t)\sin(\pi x)\sin(\pi y)$$ The jacobian then is $$F'(u)(\delta u)=\partial_t\delta u-\nabla^2\delta u$$ or in the weak form (while considering the boundary conditions) $$=\partial_t\delta u\varphi+\nabla u\nabla\varphi$$ The optimal step size should be 1, after the problem is linear.
I then compare the result for each time step with the expected results, on the one hand by calculating the $H^1$-seminorm, and on the other hand the highest value in the solution (which should be equivalent to $\exp(-2\pi^2t)$). Still, my results are not as expected:

+-------+-------------------+-----------------+
| time  | Calculated result | Expected result |
+-------+-------------------+-----------------+
|     0 | 0.336254          | 1               |
|   0.1 | 0.113069          | 0.138911        |
|   0.2 | 0.0380201         | 0.0192963       |
|   0.3 | 0.0127845         | 0.00268047      |
|   0.4 | 0.00429886        | 0.000372347     |
|   0.5 | 0.00144552        | 5.17232e-05     |
|   0.6 | 0.000486064       | 7.18493e-06     |
|   0.7 | 0.000163442       | 9.98066e-07     |
|   0.8 | 5.49583e-05       | 1.38643e-07     |
|   0.9 | 1.84801e-05       | 1.9259e-08      |
|     1 | 6.21403e-06       | 2.67529e-09     |
+-------+-------------------+-----------------+

I know that the program works for $$\partial_t u=f$$ but currently I am out of ideas about how to check why my calculated results are not the same as I would expect them. I checked the Jacobian (assumed my calculated equation above is correct), so I assume the bug is somewhere else.

Answer 1

You're using a backward Euler finite difference time stepping method. This is stable, but only first order accurate, so I suspect explains the discrepancy in the reduction factor. More explicitly, assuming a perfect spatial representation and substituting your initial condition into your difference equation gives $$(1+2\pi^2\Delta t)u^{n+1} = u^n.$$ For a timestep of 0.1, that makes each step 0.3363 of the last which (to several significant figures) is exactly what you seeing.

The only outstanding question is why your metric is giving 0.336254 rather than 1.0 for the initial condition. Based on how close that is to the number calculated above, I suspect this is actually the value for your first, rather than zeroth step.

If you want closer results, you could try using a second order (or higher) time stepping method like the midpoint method (which may or may not need surgery to your code). Note that in exchange you'll be sacrificing simplicity, stability or both.