# Tag Info

27

There is a wide variety of algorithms; Barnes Hut is a popular $\mathcal{O}(N \log N)$ method, and the Fast Multipole Method is a much more sophisticated $\mathcal{O}(N)$ alternative. Both methods make use of a tree data structure where nodes essentially only interact with their nearest neighbors at each level of the tree; you can think of splitting the ...

19

I'm not entirely familiar with what's now done for cubatures (multidimensional integration), so I'll restrict myself to quadrature formulae. There are a number of effective methods for the quadrature of oscillatory integrals. There are methods suited for finite oscillatory integrals, and there are methods for infinite oscillatory integrals. For infinite ...

18

(I have tested this approach before, and I remember it worked correctly, but I haven't tested it specifically for this question.) As far as I can tell, both $\|\mathbf{v}_1\times \mathbf{v}_2\|$ and $\mathbf{v}_1\cdot \mathbf{v}_2$ can suffer from catastrophic cancellation if they are almost parallel/perpendicular—atan2 can't give you good accuracy if ...

16

Discretization of the continuum PDEs usually commits much more error than the finite precision. I find that about 90% of the people requesting higher precision have just been lazy with problem formulation and are trying to solve a problem using poor scaling, bad discretizations, or bad continuum modeling. The remaining 10% may have justifiably ill-...

15

This question is very closely related to (and possibly a duplicate of) Is variable scaling essential when solving some PDE problems numerically?. There are still good practical reasons to nondimensionalize equations, if possible: It reduces the number of independent parameters for parametric studies (which was one of the original reasons for ...

13

Look at the fast multipole method. It is highly scalable and $O(n)$. It allows trading off between precision and cost. Here's an example where it is run at 42 Tflops on a GPU cluster.

12

Use the (IEEE standard) library function log1p, which should be present in all programming languages. The function log1p(x) returns $\log(1+x)$, and is implemented with particular attention to accuracy when $x$ is small. It is designed to solve exactly this kind of problem.

11

Besides "multidimensional vs. single-dimensional" and "finite range vs. infinite range", an important categorization for methods is "one specific type of oscillator (usually Fourier-type: $\sin(t)$, $\exp(it)$, etc, or Bessel-type: $J_0(t)$, etc.) vs. more general oscillator ($\exp(i g(t))$ or even more general oscillators $w(t)$)". At first, oscillatory ...

11

I am not an expert in machine learning, but I can outline the considerations that are relevant. The numerical calculations in machine learning are generally linear algebra -- either solving linear systems or linear least squares. For both types of problems, there are well-known backward-stable methods, so I will assume you are using a backward-stable ...

10

For $a+b\sqrt{-3}$ you can use the representation $$\begin{pmatrix}a&-3b\\b&a\end{pmatrix}$$ Addition works obviously. For multiplication, you can verify $$\begin{pmatrix}a_1&-3b_1\\b_1&a_1\end{pmatrix}\begin{pmatrix}a_2&-3b_2\\b_2&a_2\end{pmatrix} = \begin{pmatrix}a_1a_2-3b_1b_2&-3(a_1b_2+b_1a_2)\\a_1b_2+b_1a_2&a_1a_2-3b_1b_2\... 10 The numeric precision is not perfect. You get rounding errors during your computation. When working with floats, don't check if they are = 0, but check if their absolute distance to 0 is smaller than some epsilon. 10 Jean-Michel Muller, et. al., "Handbook of Floating-Point Arithmetic 2nd ed.", Birkhäuser 2018, gives the following example due to Muller, specifically constructed to deliver incorrect results with floating-point evaluation:$$ {u_{0} = 2,\\ u_{1} = -4,\\ u_{n} = 111 - \frac{1130}{u_{n-1}} + \frac{3000}{u_{n-1}u_{n-2}},\>\>\>\>n \ge 2.} $$... 9 In the 15 years that we have provided FEM software in the form of the deal.II project (http://www.dealii.org/), I don't think that we've ever had a genuine request to solve PDEs to higher accuracy than double precision. The reason is as Jed suggests in the other answer: The error one makes discretizing the PDE is much larger than the 16 digits of accuracy ... 9 Overview Good question. There is a paper entitled "Improving the accuracy of the matrix differentiation method for arbitrary collocation points" by R. Baltensperger. It's no big deal in my opinion, but it has a point (that already was known before the appearance in 2000): it stresses the importance of an accurate representation of the fact that the ... 9 I would like to know if there are any practical methods for obtaining an approximation to \mathbf{x}(t_f) (where t_f \in \mathbb{R} is some given final time) which is provably correct to N [sic] digits. That all depends on your opinion of the practicality of interval arithmetic. There are validated integrators available, such as the COSY code out of ... 9 There are aspects of modern computing systems that are inherently non-deterministic that can cause these kinds of differences. As long as the differences are very small in comparison with the required accuracy of your solutions, there probably isn't any reason to worry about this. An example of what can go wrong based on my own experience. Consider the ... 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 For the classical gravitational n-body problem, I think the following two papers do a good job at discussing the guts of the parallel implementation for the force evaluation step. Although the papers discuss a GPU implementation, they do a good job at discussing the parallelism and provide details of the algorithms: This paper by Nyland, Harris, and Prins ... 8 As an alternative source, you could also look at mesh-based Ewald-like methods. The genesis of the "particle mesh" methods (such as PPPM and smoothed particle mesh Ewald) lies in simulations of galaxies for astrophysics; the connection to charges was an unintentional side effect (that just happened to eventually overtake the original usage). More recently, ... 8 The Matlab command help eps says the following: D = EPS(X), is the positive distance from ABS(X) to the next larger in magnitude floating point number of the same precision as X. X may be either double precision or single precision. In other words, if \varepsilon_\mathsf{mach} is the relative error due to floating point, as defined in the Wikipedia ... 7 Calculate the SVD in place of the spectral decomposition. The results are the same in exact arithmetic, as your matrix is symmetric positve definite, but in finite precision arithmetic, you'll get the small eigenvalues with much more accuracy. Edit: See Demmel & Kahan, Accurate Singular Values of Bidiagonal Matrices, SIAM J. Sci. Stat. Comput. 11 (1990)... 7 For all non-trivial problems (i.e., for those where performance matters) almost all of the memory you have will be in the matrix, and relatively little in vectors. For example, for 3d Taylor-Hood elements for the Stokes equation, you have a few hundred elements per row in the matrix, and this vastly outweighs the amount of memory needed for vectors. We have ... 7 It is the difference in scales between terms in your equations that tend to cause numerical difficulties. You may work in any units you like as long as you are consistent. My approach has been to always consistently non-dimensionalize my equations in order to reduce the number of parameters to the minimum required, but this is only for my convenience. 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 [EDIT] An alternate view: 64-bit floating numbers represent a discrete set S. For a function f to be exactly invertible, it should be a bijection from S to S. Suppose we are interested in a fast growing function like \exp. At one point, \exp x > x. If M is the maximum element from S, then f(M)\ge M. The strict version f(M)> M is ... 7 The efficient answer to this question is, not too surprisingly, in another note by Velvel Kahan:$$\alpha=2\arctan\left(\left\|\frac{\mathbf v_1}{\|\mathbf v_1\|}+\frac{\mathbf v_2}{\|\mathbf v_2\|}\right\|,\left\|\frac{\mathbf v_1}{\|\mathbf v_1\|}-\frac{\mathbf v_2}{\|\mathbf v_2\|}\right\|\right) where I use $\arctan(x,y)$ as the angle made by $(x,y)$ ...

7

Added after my initial answer: It appears to me now that the author of the referenced paper is giving condition numbers (apparently 2-norm condition numbers but possibly infinity-norm condition numbers) in the table while giving maximum absolute errors rather than norm relative errors or maximum elementwise relative errors (these are all different ...

6

This may be a bit outdated, but Hairer and Wanner's book recommends their own radau5 code. In my own experience, this code is both extremely robust and efficient. It is also a rather straight-forward Runge-Kutta scheme, and thus not too difficult to implement. There's a Matlab implementation of radau5 by Christian Engstler somewhere out there.

6

There are many good open-source implementations of root-finding methods out there already. One example is boost, whose implementation of Newton-Raphson and related methods you can find here. If you read its source code, you will be able to see what issues the authors wanted to address, and how they dealt with issues like convergence, user-specified tolerance ...

6

There are three issues that are likely to cause such problems in pseudospectral methods: Gibbs oscillations Aliasing Time step too large In any case you likely develop oscillations in the solution until some point ends up with a negative density, resulting in a NaN when computing the pressure or sound speed or some other term. The solution to 3 is obvious, ...

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