0
$\begingroup$

In MATLAB's ODE suit, there seem to be two basic methods of controlling the Local Truncation Error (LTE) of the ODE which the user can choose from, namely:

  1. The absolute error control (default), |e(i)| <= max(RelTol*abs(y(i)),AbsTol(i))
  2. The norm error control (setting option normcontrol, on), where norm(e(i)) <= max(RelTol*norm(y(i)),AbsTol(i))

By alternating between the two for the same problem, in the case of using the ode23tb for example, the difference seems to be noticable as seen in the figure below. Furthermore, in the case of the norm error control the algorithm is significantly faster as well.

enter image description here

Thus my question is threefold:

  1. The math may be apparent by its formulation, but is there an intuitive physical interpretation of the norm error control method?
  2. What are the advantages and disadvantages of using the norm error control option as compared to the absolute error?
  3. In which context could the one method be more accurate than the other?
$\endgroup$
1
  • $\begingroup$ The documentation is written in a confusing way, and doesn't fully specify what is done. However, it seems that the default uses the maximum norm while "norm control" uses the $\ell_2$ or Euclidean norm. Which is preferable will depend on what the solution actually represents. Which one is more accurate would depend on how the $\ell_2$ norm is defined (i.e., whether there is some normalization with respect to the dimension of the problem). But you can adjust the accuracy anyway by changing the tolerances. $\endgroup$ – David Ketcheson Dec 6 '20 at 11:31
1
$\begingroup$

As far as I understand it, if you do not use the NormControl, the time step is adapted so that the maximum error on any of the solution components is below the tolerance threshold.

If, on the other hand, you use the NormControl option, then the time step is adapted so that the overall error norm is lower than the tolerance threshold. This is slightly less stringent as the biggest soltion component error is "miwed" with the other errors. Let's take a worst-case scenario: you have $N-1$ variables that follow an ODE of the type $y' = a$ with $a$ a constant. Runge-Kutta methods integrate this exactly so the error on these components is zero to machine precision. Now if you had one solution component (the $N$-th variable) which follows a nonlinear ODE, the integration error estimate $e_n$ on this variable will be nonzero. Without the NormControl option, the time step will be adjusted roughly as follows: $$ \Delta t \approx \left( \dfrac{|e_N|}{atol + rtol |y_N|} \right)^{\frac{1}{q+1}}$$ with $q$ the order of the integration error estimate, which is most likely equal to $p-1$, with $p$ the order of the method.

With the NormControl option set, it will be adjusted as: $$ \Delta t \approx \left( \dfrac{ \sqrt{\sum\limits_{i=1}^{N} |e_i|^2}}{atol + rtol \sqrt{\sum\limits_{i=1}^{N} |y_i|^2}} \right)^{\frac{1}{q+1}}$$

which, in our worst-case scenario, would be: $$ \Delta t \approx \left( \dfrac{ |e_n|}{atol + rtol \sqrt{\sum\limits_{i=1}^{N} |y_i|^2}} \right)^{\frac{1}{q+1}}$$ If we compare this to the very first equation of this answer, we see that the numerator is the same, but the denominator is larger, as $||y|| \geq |y_n|$. So the integration error will be perceived as lower with NormControl, hence the time step will be larger and the simulation quicker. But this also means that the solution quality is slightly worse.

So in your figure, I think the solution without NormControl is the better one in terms of solution exactitude. You could check that by lowering the tolerance down to very fine levels (say 1e-12) to see what the "exact" solution is. Also, you can compare how the time step evolves between the two cases to get more insight on this behaviour.

I am not quite sure, but I think the error control based on the norm of the error may be more gentle with the time step variations, whereas the solution based on the maximum error component may have bigger "jumps" in the error estimate as the solution evolves.

Also, I usually use the norm $||y|| = ||y||_{Matlab} / \sqrt{N}$, i.e. a root mean square error. This way, for discretized PDEs, the error does not "fictiously" grow as the number of mesh points increases. However, single large error components may not be represented well as they will be averaged with the lower components.

The norm in Matlab does not seem to have the factor $1/\sqrt{N}$. Therefore, this should mean that, if your system is sufficiently large, the NormControl option will actually be more stringent on the time step. Indeed, say we duplicate a scalar ODE $N$ times, then the numerator of the above time step formula will be: $\sqrt{\sum\limits_{i=1}^{N} |e_i|^2} = \sqrt{N |e_1|^2} = \sqrt{N} |e_1|$ which will be higher than $max(|e|)$. Thus I am guessing that your test case had a number of variables $N$ reasonably low.

In any case, the choice of using or not using this NormControl option should not affect wildly your result, otherwise this means that your integration tolerances are too large.

$\endgroup$
1
  • 1
    $\begingroup$ It should be $\frac1{q+1}$ in the exponents for the next step size, where $q,q+1$ are the orders of the embedded method and the method step is the one of the order $q+1$ method. $\endgroup$ – Lutz Lehmann Dec 6 '20 at 22:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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