# Debugging Newton-method used in a CG-approach

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.

• Are you discretising for time in a finite element sense? – origimbo Feb 1 '19 at 12:50
• I added my time discretization in the question – arc_lupus Feb 1 '19 at 12:53

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.