Asymptotic error of forward Euler

Question

I'm trying to understand the asymptotic error behaviour of forward Euler (finite difference method), as timesteps are decreased (refined), so I feel trust in the method of manufactured solutions (MMS).

But looking at a simple model of a ball accelerating, I noticed you can make the error asymptote to whatever you want...

If the ball begins with zero velocity, then forward euler estimates that it does not move in the first timestep. So the error is the same as the analytically determined position:

• Constant acceleration gives velocity $v(t) = t$ and position $p(t)=\frac{1}{2}t^2$. So that is the error, which is order $O( (\Delta t)^2 )$, for the first timestep to $t=\Delta t$. So a plot of error vs timestep, with timestep refinement, will reveal an asymptote of that order.

• An increasing acceleration, say $v(t) = t^2$, gives $p(t)=\frac{1}{3}t^3$. The error is now order $O( (\Delta t)^3 )$, which it again asymptotes to with timestep refinement.

This is troubling, because MMS relies on error vs timestep having the same asymptotic behaviour of the estimation method (in this case, forward Euler, which is a first order method). In this simple example, it doesn't seem to work.

MMS is a well-established technique, so I must be doing something wrong or misunderstanding something...

1. Is it that you should take many timesteps (accumulating error), not just the first one?

2. Is it that this example has only one variable $p$, but fluid simulations have many (a grid of many cells, with many variables in each), and in a complicated scenario, with typical fluid equations, it kind of all evens out?

maybe relevant: Numerically determining convergence order of Euler's method

---

To summarize my understanding (seeking feedbac!):

The error term expected to dominate evaluates to zero

Forward Euler being first order is confusing to me, because for estimating the function, there's local truncation error (LTE) $O(h^2)$ and global trunction error $O(h)$, and for estimating the first derivative (used in forward Euler/forward differences) the order of accuracy is $O(h)$. Here, I'll stick with LTE because it's applicable to the single timestep used here and (bonus!) I kinda understands it.

In the Taylors series expansion of a function, the first two terms can be seen as forward Euler, $p(t+h) \approx p(t)+hp'(t)$, with the remaining terms being trunction error (just two error terms shown here):

$ p(t+h) = p(t)+hp\prime(t) + \frac{h^2}{2} p\prime\prime(t) + \frac{h^3}{3} p\prime\prime\prime(t) + \ldots $

The LTE order of the error is $O(h^2)$, which will be greater than the higher order terms, because $O(h^2)<O(h)$ for $h \rightarrow 0$ (similar to why $x^2<x$ for $x<1$).

Applying to the two examples above:

• $p(t)=\frac{1}{2}t^2$, with derivatives $t, 1, 0$, at $t=0$ are $0,1,0$: $p(t+h) = 0 + h.(0) + \frac{h^2}{2}.(1) + \frac{h^3}{6}.(0)$ Here, the truncation error is $\frac{h^2}{2}$, which is order $O(h^2)$ - exactly as expected.

• $p(t)=\frac{1}{3}t^3$, with derivatives $t^2, 2t, 2$, at $t=0$ are $0,0,2$: $p(t+h) = 0 + h.(0) + \frac{h^2}{2}.(0) + \frac{h^3}{6}.(2)$ Unfortunately, the truncation error is $\frac{h^3}{3}$, which is order $O(h^3)$ - NOT as expected. The problem is the error term expected to dominate ($\frac{h^2}{2}.(2t)$) evaluates to zero (at $t=0$).

This is one danger warned of in the answer by @cpraveen: zero derivatives in the truncation error. I read this at first as meaning the zero function (like the comment by @wolfgangbangerth, $p\prime(t)=0$), but it happens here because the derivative expected to dominate evaluates to zero. A derivative evaluating to zero is also a zero!

Also (if I understand correctly) non-linearity in the answer would be a problem because its higher derivatives are zero (similar again for non-trivial in the answer). An advantage of $e$ and $\sin$ in manufactured solutions is they never differentiate to zero... though all $\sin$ derivatives evaluate to zero twice per cycle, so that might cause the same problem, if that special input values come up too often (e.g. wavelength aligned with grid). This fits with variation in the solution.

Avoiding it in the examples

Simply evaluating at $t=1$ (instead of $t=0$) fixes it, because the error term expected to dominate no longer evaluates to zero.

• $p(t)=\frac{1}{2}t^2$, with derivatives $t, 1, 0$, at $t=1$ are $1,1,0$: $p(t+h) = 1/2 + h.1 + \frac{h^2}{2}.1 + \frac{h^3}{6}.0$ The truncation error is still order $O(h^2)$ as before (good!)

• $p(t)=\frac{1}{3}t^3$, with derivatives $t^2, 2t, 2$, at $t=1$, are $1, 2, 2$ $p(t+h) = 0 + h.1 + \frac{h^2}{2}.2 + \frac{h^3}{6}.(2)$. The truncation error now has both terms non-zero, so the order is $O(h^2)$, which is what it should be. (fixed!)

tl;dr Conclusion

It seems the cause is just these special cases, of zero error terms. One way to avoid this is to use $e^t$ in the manufactured solution, because its derivatives never evaluate to zero. Another is to use $\sin(t)$, though its derivatives do evaluate to zero sometimes.

In a proper fluid simulation, the interactions between equations could well produce zero error terms from time to time, no matter what you do. It seems the only way to guard against this is to have varied manufactured solutions, so that zero error terms only occur rarely in the grid, and are overwhelmed by the error terms of all the other grid cells. (But careful of systemic zeros, like from cell-aligned $\sin$ periods). (Using many steps instead of one timestep would also help, for the same reason.)

We must remember that error is supposed to be random, so that this probabilistic averaging out is perfectly OK!

BTW It reminds me of the probablistic Schwartz-Zippel lemma: (oversummarizing) you can test for the zero polynomial by evaluating with random arguments. But it could be zero by chance; that is, you hit a root - a special case. The lemma shows the probability of this happening, so if you try several (unrelated) random arguments, you can make the probability as small as you like.

I guess a similar analysis might have been done for MMS - but it seems so much more difficult because of the complexity of the governing equations, over a grid, over time.

Have I got all that about right?

---

I did some numerical tests, which comfirmer the above. Surprisingly (to me), for $p(t)=\frac{1}{3}t^3$, at $t=1$, I didn't get exactly $2$ for the power law, but around $2\pm .02$. I only played around a little, but seemed to be that, even with much smaller refinements in $h$ of $0.001$.

I have heard that that level of accuracy is pretty good for MMS.

I also tried $p(t)=\sin t$ and $p(t) = e^t$, and got that $2\pm.02$ power law again.

And I confirmed that $p(t)=\sin t$ doesn't give the expected err vs timestep asymptote of $2$, for $t=0$, because the second derivative is $0$ there. But it works for $t=1$.

Further note: the error term expected to domimate being zero is a particular case, of that term being smaller than following terms. So, for $\sin t$ this would also happen in s small region around $t=0$ not just exactly at zero. How small a region depends on how small $h$ is, because the higher powers of $h$ in the following terms makes them smaller.

I think we can see all this as a natural probabilistic distribution of random errors. $O(h^2)$ is approximate, because there are still other error terms. Because special cases are rare, if they sometimes do occur, it doesn't really matter - provided there is adequate sampling of the overall distribution, to even them out. I guess that explains the $2\pm0.02$, and it's not just this special case, but general random variation in the errors.

It seems my strategy, of starting with the simplest possible case of a sample of one from the distribution, was bound to have problems! EDIT but much easier to analyse.

But we could attempt to minimize these special cases by choosing non-zero regions of $\sin$ - especially for time $t$, which would be used throughout a simulation grid. It might also be helpful to similarly choose a non-zero region of $\sin$ over the entire grid (i.e so that $\sin x$ is zero outside the grid, e.g. $x=0$ before the grid starts, and $x=\pi$ after the grid ends).

---

I tried detecting mistakes with this, since that's the purpose of MMS.

For forward Euler, the derivative needs to be non-zero, or that term won't be exercised, and any mistakes in it won't be detected. For $p(t+h) \approx p(t)+hp'(t)$, we don't want $p\prime(t)=0$. This happens e.g. for $p(t) = \sin t, t=\pi/2$.

Using $p(t) = e^t$ avoids this problem, but won't detect using the wrong order of derivative (because they're all the same). Can fix this with e.g. $p(t)=e^{2t}$.

I also found gnuplot fit had problems with $p(t) = e^{t^2}$, when $t>4$, with a power law $ax^b$, but was better if converting to logs before fitting (and using a linear law - how you'd do it by hand).

Answer 1

When we say that Euler method is first order accurate, it means that for a class of ode with sufficiently smooth solutions, the error will be at most $O(h)$ and there is at least one ode for which it cannot be smaller than $O(h)$. So it can happen that for some ode, the error will be smaller than $O(h)$.

The example you give is somewhat special. But that also happens at the first time step only.

Hence when using MMS, it is important to take non-trivial solutions which have a lot of variation and non-linearity and have all derivatives non zero that may appear in the truncation error. MMS is a way to verify theoretical convergence rates and code correctness, and it is not a proof of the convergence rate.

Answer 2

MMS for first order ODE

The method of manufactured solutions takes some ODE $$\dot x(t)=f(t,x(t))+r(t)$$ (where $\dot x=\frac{dx}{dt}$) and adapts the perturbation or control term $r$ in such a way that some previously selected function $p$ is a solution. This means setting $$r(t)=\dot p(t)-f(t,p(t))$$ and $x(t_0)=p(t_0)$. For instance, for $$\dot x+\sin(x)=r~~\text{ and }~~p(t)=\sin(2t)$$ one would have to set $$x(0)=0~~\text{ and }~~r(t)=2\cos(2t)+\sin(\sin(2t)).$$

---

The error of the Euler method has the form $O(t\cdot Δt)$ for small $t$ or a little more globally $O((e^{Lt}-1)⋅Δt)$ where $L$ is a Lipschitz constant close to the initial point (see Grönwall lemma and that the numerical solution can be considered an $O(Δt)$ perturbation of the original ODE).

Using the above example over the interval $[0,3\pi]$ and dividing the (absolute) error of the Euler method by $Δt$ to get comparable results for the different step sizes result in the following plot The zero crossings of the error function are caused by the oscillating nature of the (signed) error. However, apart from these points, one can see the convergence to a constant (relative to $Δt$) coefficient of the leading error term $C(t)Δt$. And the global error is related to the worst or maximimum local error, however that is measured.

What is wrong in the approach in the question

If you do an asymptotic error analysis, you do not observe the error in the first step of the varying $Δt$, but the error at some fixed time $t$ that requires $n=n(Δt)=(t-t_0)/Δt$ integration steps to reach.

The ODE you seem to solve as test case is simply a quadrature problem $\dot x(t)=r(t)$, where then consequently for MMS $r(t)=\dot p(t)$ for some given test function $p(t)$. Then for all $Δt$ chosen such that exactly $t=t_0+nΔt$ for some integer $n$ the Euler method is simply the Riemann sum with left sample points, where one can find a more concrete form of the leading error term (see also the Euler-Maclaurin summation formula) \begin{align} x_n&=p(t_0)+\sum_{k=0}^{n-1}\dot p(t_0+kΔt)Δt \\ &=\sum_{k=0}^{n-1}[p(t_0+(k+1)Δt)-p(t_0+kΔt)+\frac12\ddot p(t_0+(k+\theta_k)Δt)Δt^2] \\ &=p(t_0+nΔt)+\frac12\sum_{k=0}^{n-1}\ddot p(t_0+(k+\theta_k)Δt)Δt^2 \\ &=p(t_0+nΔt)+\frac12[\dot p(t_0+nΔt)-\dot p(t_0)]Δt +O(Δt^2) \end{align} That last expression tells us that the error is always $O(Δt)$, with the leading error term vanishing only for $\dot p(t)=\dot p(t_0)$. This can be true for all $t$ only if $p(t)$ is a linear function where the Euler method integrates the constant derivative exactly. In testing the method, one should avoid the roots of $\dot p(t)-\dot p(t_0)=0$ where the errors approximately cancel systematically, as that may be specific only for quadrature problems and not for more general, proper ODE.