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


25

Take \begin{align} 1-\sqrt{ 1-x^2} &= (1-\sqrt{ 1-x^2})\frac{1+\sqrt{ 1-x^2}}{1+\sqrt{ 1-x^2}}\\ &= \frac{x^2}{1+\sqrt{ 1-x^2}} \end{align} So \begin{align} y = x\sqrt{\frac{1}{2+2\sqrt{1-x^2}}} \end{align}


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


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


13

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

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


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


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 λ&...


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

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


9

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

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

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

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


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

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

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

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


4

I am not sure this is possible with the Python libraries since they are using Fortran under the hood and that can't be easily recompiled, but the Julia DifferentialEquations.jl JIT compile specializes the solvers based on the number types that you give it. Here's a demonstration of some weird types like rational numbers, MPFR BigFloats, and ArbFloats (based ...


4

Stability does not necessarily imply accuracy. I'll demonstrate this with a simple scalar ODE $y' = \lambda y$ (known as Dahlquist's test equation). This simple ODE is generally interpreted as a linearisation of a generic ODE $y'=f(y)$ around a given solution point, i.e. it reproduces, at least locally, the behaviour of the more complex ODE you might be ...


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

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

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


3

This is what Griewank et al. call "Piecewise linearization in secant mode", see for instance https://opus4.kobv.de/opus4-zib/files/6164/newton_secant_approx_paper.pdf. The aim of that research was to capture the kinks of absolute value operations with the same precision a tangent or a secant captures the local behavior of a smooth function, with an ...


3

The choice of finite-difference scheme depends on several factors, such as the smoothness of your data, how uniformly-spaced the data actually is, etc. You may also want to consider just how accurate your velocity estimate actually needs to be. For example, if there are large error bars in the trajectory then it probably makes little sense to use a high-...


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

When solving a linear system $Ax=b$, you suffer from two sources of error: Round-off error The ill-conditionedness of the matrix. The condition number of a computational representation of a matrix is, for all practical purposes equal to the condition number of the exact matrix (which you can't represent computationally), so it does not matter whether you ...


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