3
$\begingroup$

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).

$\endgroup$
  • 1
    $\begingroup$ yes you summarised it well. Such things happen rarely when you go to real, complicated problems. If you take non-linear, non-polynomial solutions, which have variations in all independent variables, and enough derivatives including mixed ones are non-trivial, then you should be ok to use mms. $\endgroup$ – cpraveen Nov 23 '18 at 6:15
  • $\begingroup$ @cpraveen Thank you! I just noticed your comment now, which you posted 37 min ago. I probably added a fair bit since you saw it. $\endgroup$ – hyperpallium Nov 23 '18 at 6:54
4
$\begingroup$

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.

$\endgroup$
  • 3
    $\begingroup$ In particular, all useful methods will have a zero error if the ODE happens to be $\dot x(t)=0$. This does not mean that they will have a zero error for all ODEs, of course. $\endgroup$ – Wolfgang Bangerth Nov 20 '18 at 14:03
  • $\begingroup$ Thanks. Is that sort of my guess number 2 ("it all evens out")? (I edited numbering in) I've now read more, and examples seem to use a fixed time $t$ to compare the different timesteps (always more than a single timestep, and increasing as $\Delta t$ decreases), so I was thinking my guess number 1 was also true? (Or that just for easy comparison?) $\endgroup$ – hyperpallium Nov 21 '18 at 7:37
  • $\begingroup$ What would you recommend for building trust in MMS? My goal isn't to prove convergence rates, but to be able to check my codes with MMS. For that, I need to have confidence in MMS... how could I apply your advice (variation, non-linearity and derivatives) to this simplest example of forward Euler, so that MMS works? $\endgroup$ – hyperpallium Nov 21 '18 at 7:44
  • 1
    $\begingroup$ When you have both spatial and time errors, and you want to test accuracy of time scheme, then doing just one step is still useful. If spatial error are dominant, then only reducing time step may not give enough reduction in error to quantify time integration errors. But you have to choose somewhat non-trivial solutions in this case, and do many of them to gain confidence. Maybe try non-polynomial solutions. $\endgroup$ – cpraveen Nov 21 '18 at 7:56
1
$\begingroup$

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 error 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.

$\endgroup$
  • $\begingroup$ Thanks for confirming more than one timestep is the way to go. I googled and read up on your answer, but it's over my head for now. Some quick questions: What does $\dot x=r=\dot p$ mean? (I'd guess: the derivatives of functions $x$ and $p$ equal constant $r$, but I bet that's not right) I think that avoiding the roots of $\dot p(t)-\dot p(t_0)=0$ would be really difficult to do, as it's hard to even know what they are? And... wouldn't it also be a problem if $\dot p(t)-\dot p(t_0)$ was much smaller than the $O(Δt^2)$? (i.e. not just for equalling zero) $\endgroup$ – hyperpallium Nov 26 '18 at 8:41
  • 1
    $\begingroup$ Yes, the sample time should stay constant, not the number of time steps. -- This is how MMS works, you seem to solve $\dot x=\frac{dx}{dt}=r$ where you select $r$ so that $x=p$ is a solution, thus $r=\dot p$. Meaning that no, $r$ is not a constant, just the right side. -- As the test function $p$ is selected by you, you also have full control over the derivative values. -- No, $\dot p(t)-\dot p(t_0)$ is a constant, $O(1)$, while $Δt$ is variable for this experiment. $\endgroup$ – LutzL Nov 26 '18 at 9:21
  • $\begingroup$ Thanks. Going to my 2nd question of 3 in my comment: $t_0$ is the start time, and $t$ is the end time - both are constants. The roots of $\dot p(t)-\dot p(t_0)$ are the zeros, i.e. where $\dot p(t)-\dot p(t_0)=0$. To find those roots, you'd substitute the specific $t_0$ and $t$ into the derivative, and subtract. I suppose that might or might not be difficult! For the first example in my Question, of $p(t)=\frac{1}{2}t^2$, to have a root .you'd have to choose $t_0=-t$, e.g. $t_0=-1, t=1$. But in general, there could be many roots... $\endgroup$ – hyperpallium Nov 27 '18 at 8:41
  • 1
    $\begingroup$ Yes, the leading error term is $(\dot p(t)-\dot p(t_0))Δt$, and for $t$ where this coefficient is zero you get that the error at that point, and only at that point, is $O(Δt^2)$. $\endgroup$ – LutzL Nov 27 '18 at 9:09
  • 1
    $\begingroup$ Then you would have to make some restrictions on $Δt$ to not be arbitrarily small. And also compute the coefficients to some accuracy. An error $0.1Δt+ Δt^2+..$ will look quadratic for $Δt>0.1$ but is linear for smaller $Δt$. In the general asymptotic analysis $Δt$ is to be considered arbitrarily or infinitely small, and thus it is only important if coefficients are zero or not. $\endgroup$ – LutzL Nov 27 '18 at 9:24

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.