35

Since I just finished optimizing a lot of them in a software, DifferentialEquations.jl, I decided to just lay out a comparison of the main Order 4/5 methods. The Fehlberg method was left out because it's commonly known to be less efficient than the DP5 method. Backstories Dormand-Prince 4/5 The Dormand-Prince method was developed to be accurate as a 4/5 ...


14

Here is R1, as computed in MATLAB: 1.0e+07 * -7.382605957465515 -9.599867106092937 -2.830412177259742 -0.000000000002830 -0.000000000002830 -1.230434326244253 -1.599977851015490 -0.471735362876624 -0.000000000000472 -0.000000000000472 3.691302978732758 4.799933553046468 1.415206088629871 0.000000000001415 0.000000000001415 -5....


14

First, see Mark L. Stone's answers, which is completely correct. Second, realize that this is the reason why people told you to use relative errors in your numerical analysis class. :) Third, the real question here is why the results do not coincide exactly, since both languages call some BLAS library functions for their computations. There are several very ...


12

Here is a quick example which is very practical related to memory usage in PDEs. When one discretizes the Laplace operator, $\Delta u$, for example, in the Heat Equation $$ u_t = \Delta u + f(t,u) .$$ To solve it numerically, one ends up with sparse matrices $A$, and a method of lines discretization then solve $$ u_t = Au + f(t,u) $$ The canonical 1D ...


12

I was able to reproduce the behavior reported in the question, and traced the observed inaccuracies to the following line: return y*sin(pi<Real>()*x)/pi<Real>(); The explicit multiplication with a floating-point approximation of π introduces a small error into the argument to sin, which comprises the representational error in the constant and ...


11

Normally there are some principal reasons to prefer solve a linear system respect to use the inverse. Briefly: problem with the conditional number (@GoHokies comment) problem in the sparse case (@ChrisRackauckas answer) efficiency (@Kirill comment) Anyway, as @ChristianClason remarked in comments, can be some cases where the use of the inverse is a good ...


11

$p \log p$ won't suffer from precision loss anywhere in $[0,1]$, and won't suffer from exponent overflow near $0$ either. Thus, the fast, accurate way is p ? p * -log(p) : 0


9

The integral in question is also known as the Boys function, after the British chemist Samuel Francis Boys who introduced its use in the early 1950s. A few years ago, I needed to compute this function in double precision, as fast as possible but accurately. I managed to achieve a relative error on the order of $10^{-15}$ across the entire input domain. It ...


9

You can sometimes prove such results (or get counterexamples) using an SMT solver such as Z3 that supports floating point arithmetic. Here is a proof of a version of your theorem that says $|((x+y)-y)-x| \leq 2^{-23}|x|$ when $x>y>1$ and $x+y\neq\infty_{32}$ in 32-bit floating point arithmetic: λ> import Data.SBV λ> :set -XScopedTypeVariables λ&...


8

It's fairly easy to evaluate, to do this expand the logs in Taylor series in $x=(k+1)^{-2}$: $$ \log_2(1-x) = \frac{-1}{\log 2}\sum_{m\geq1}\frac{x^{m}}{m}$$ $$ \log_2(-\log_2(1-x)) = \frac{\log x}{\log 2} - \frac{\log\log 2}{\log 2} + \sum_{n\geq 1}a_n x^n, $$ where $a_n$ are Taylor series coefficients of the l.h.s. after the log-singularity is subtracted. ...


7

At least in MATLAB, I believe abs(z) is implemented as sqrt(z*z'). The extra square-root and squaring operation reduces numerical precision. >> z = randn + randn * i z = 0.5377 + 1.8339i >> abs(z)^2 - z*z' ans = 4.4409e-16 >> abs(z)^2 - sqrt(z*z')^2 ans = 0


7

The inverse Langevin function $\mathcal{L}^{-1}(x)$ is an odd function. Therefore one needs to consider only approximation on the interval $[0, 1]$; the negative half-plane is treated via symmetry about the origin. Because of the singularity at unity, a polynomial approximation is not a fruitful direction of the investigation, although some authors have ...


6

I think you probably already know all this, and maybe it's just the wording that is confusing. I'm going to rephrase what they're doing to make it more explicit. It's exactly the same calculation of truncation errors that you are already familiar with. When approximating an operator $\mathcal{L}$ with an FD operator $L\approx \mathcal{L}$, the truncation ...


6

For a 6x6 grid, those are about the error differences I would expect from two different methods. You have to realize that a 6x6 grid is a very coarse grid, even for a simple problem like yours. As long as you see the two solutions converge towards each other as you refine your grid, there is likely no implementation error. Finite-difference has no general ...


6

Use equation (2), equation (1) is wrong. Technically the $L^2$ norm (upper case "L") is an integral norm of a function defined as $$ \left|\left|f(x)\right|\right|_2 = \sqrt{ \int_\Omega |f(x)|^2 dx}. $$ I'm sure you meant $l^2$. What you typically want to compute in the context of convergence of numerical methods is the finite dimensional analog of the $L^...


6

If you want to increase the accuracy of a finite difference scheme, you can always try increasing the degree of your stencil. On equidistant points, though, this can lead to numerical instabilities. To avoid these problems and still get high accuracy, I would suggest using Spectral Methods. If your problem has fixed poles, you can try to get around them by ...


5

You might take a look at Numerical Methods for Special Functions by Amparo Gil, Javier Segura, and Nico M. Temme.


5

The purpose of automatic error estimation and step size control is to free you from the problem of determining manually what a sufficiently small step size is. So your question is a bit like asking "somebody gave me this automatic transmission car; how can I tell what gear I'm in?" The point is that you shouldn't need to know. Of course, if the ...


5

According to [1]: "However, the IEEE-754 standard specifies nothing for elementary functions" and "Indeed, the mathematical libraries (libm) provided by operating systems do not guarantee correct rounding.". As such, according to [1], the answer is that the IEEE-754 standard does not require the exp2 function to be correctly rounded, it only recommends it. [...


5

Short version In scientific computing, the notion of relative error is way more popular than accuracy to $N$ digits. Whenever we present the results, we usually plot the obtained (scaled) data and relative/absolute error in multiple different ways. Reporting the number of correct obtained digits (in addition to the relative error) is usually limited to the ...


5

What you are looking for is called "bootstrapping". It is a common problem of all multistep ODE integrators and is discussed in many books on the topic. Among your options are to use a lower-order method with smaller time step, or to use a one-step method of higher order for the first few steps (e.g., a Runge-Kutta method).


4

I'd take a look at Abramowicz & Stegun's book, or the newer revision that NIST has published a couple of years ago and that's available online I believe. They also discuss ways to implement things in a stable way.


4

Yes. You can compute a running error bound, i.e, a number $\mu$ such that the difference between the exact value of $y = p(x)$ and the computed value satisfies $\hat{y}$ satisfies $$|y - \hat{y}| \leq \mu u.$$ Here $u$ is the unit roundoff. You can trust the sign computed sign of $y$, when $|\hat{y}| > \mu u$. Let $p(x) = \sum_{j=0}^n a_j x^j$, then ...


3

Your initial method for computing the norms was incorrect (see my answer to your related question) The norm was incorrect by a factor of $1/\sqrt{N}$, leading to underestimating the error as $N$ increases and giving the appearance of a higher-than-expected convergence rate. I would expect that both L1 and L2 norms give somewhere between orders 1 and 2 for ...


3

Compensated Horner method (http://www-pequan.lip6.fr/~jmc/polycopies/Compensation-horner.pdf) has an error bound of the form $$ |\mathrm{comphorner}(p, x) - p(x)| \leq u|p(x)| + \gamma_{2n}^2\tilde p(x), \qquad \tilde p(x) = \sum_i |a_i||x|^i, \quad \gamma_n = \frac{n u}{1-n u}, $$ where $u$ is the unit roundoff, so it will give the correct sign so long as $|...


3

If you have theoretical expectations for some observables of simulations, I see two general ways of dealing with them: You exploit them to get more accurate results, e.g., you make your algorithm use the fact that $F_1=F_2$. In this case, the value of your observables is not interesting since it they match your theoretical result by construction. You use ...


2

So I went ahead and implemented a code in Matlab that can solve this problem using a spectral approach, utilizing a simple polynomial basis. Using a simple polynomial basis of order 20 resulted in the following: One thing I did do that you might not have was define the basis to exist on a domain that might have better numerical stability (especially as x ...


2

This depends on the problem you are solving. Usually the elliptic and hyperbolic cases are handled separately because the spectral properties of their operators are very different, and so you have papers/books which answer this question for those specific cases. Generally speaking the expectation is that if you use p-th order polynomials, then DGFEM/FEM ...


2

You could use interval arithmetic for this purpose.


Only top voted, non community-wiki answers of a minimum length are eligible