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

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

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

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

$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

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


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


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


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


5

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


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


5

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

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

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


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