33

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

The best solution that I know of is to program the symbolic expressions in Mathematica, Maple, or SymPy; all of the links go directly to the code generation documentation. All of the programs above can generate code in C or Fortran. None of the programs above mentions accuracy in IEEE 754 arithmetic; in general, it would be difficult to anticipate all ...


14

If you want an idea of just how far we are away from such a software package, please look at the 2001 LAPACK working note on computing Givens rotations reliably and efficiently. I would expect most non-specialists (and many specialists!) in numerical analysis to be surprised at just how much analysis went into solving such an ostensibly simple problem: ...


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

Code generation and precompilation of mathematical expressions is becoming more popular. While symbolic packages like SymPy, Mathematica, and Maple may include code generation I'm not confident that any of them also think hard about numerics. There are a couple of other projects one could look into that are interested both in symbolics and in numerics. ...


11

In practice, most people stick to relatively low orders, usually first or second order. This view is often challenged by more theoretical researchers that believe in more accurate answers . The rate of convergence for simple smooth problems is well documented, for example see Bill Mitchell's comparison of hp adaptivity. While for theoretical works it is ...


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


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


10

This can be done by integration by parts: $$ \int^\infty_1 x e^{-ax} = \frac{-1}{a} x e^{-ax}\mid^\infty_1 - \frac{-1}{a} \int^\infty_1 e^{-ax} = \frac{e^{-a}}{a} + \frac{e^{-a}}{a^2} = \frac{a+1}{a^2} e^{-a} $$ and continuing on by induction $$ \int^\infty_1 x^k e^{-ax} = \frac{-1}{a} x^k e^{-ax}\mid^\infty_1 - \frac{-k}{a} \int^\infty_1 x^{k-1} e^{-...


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


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

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

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

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

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

As far as finite element solid mechanics simulations are concerned, you can't use less than 8 quadrature points without using stabilization forces. In case of incompressible material (your case), the best solution for accuracy purpose is to use mixed formulation. You can refer to the book by Simo and Hughes : http://books.google.fr/books/about/...


4

Not in general, I can safely say the implementor of the code generator in SymPy didn't even try =P. Paolo Bientinesi developed a method for generating stability proofs of linear algebra algorithms, which are generated using Robert van de Geijn's FLAME notation. See this paper, or a longer, working note version.


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

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


3

Guidelines: High-order methods for problem where solution is expected to be smooth and otherwise low-order methods and/or methods which can handle discontinuities in solutions. In cases where high-order methods can be exploited there can be significant saving in computational effort measured in terms CPU time as a result of high convergence rate. For ...


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