Kirill
  • Member for 10 years
  • Last seen more than a week ago
How to avoid catastrophic cancellation in python function?
Accepted answer
32 votes

This is indeed called catastrophic cancellation. In fact, this particular case is very easy: rewrite the function using the equivalent, numerically stable expression $$ \frac{t}{1+\sqrt{1-t^2}}. $$ ...

View answer
Numerically stable way of computing angles between vectors
Accepted answer
22 votes

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

View answer
About faster approximation of log(x)
Accepted answer
21 votes

This is not really an authoritative answer, more a list of issues I think you should consider, and I haven't tested your code. 0. How did you test your code for correctness and speed? Both are ...

View answer
Robust computation of the mean of two numbers in floating-point?
18 votes

I think Higham's Accuracy and Stability of Numerical Algorithms addresses how one can analyze these types of problems. See Chapter 2, especially exercise 2.8. In this answer I'd like to point out ...

View answer
Super C++ optimization of matrix multiplication with Armadillo
Accepted answer
14 votes

In fundamental C++, I find the problem here is that C++ will allocate a new object of cx_mat to store evolutionMatrix*stateMatrix, and then copy the new object to stateMatrix with operator=(). I ...

View answer
Eigenvalues of Small Matrices
Accepted answer
13 votes

The first thing to note is that the correspondence between finding roots of a polynomial (any polynomial) and finding the eigenvalues of an arbitrary matrix is really direct, and it's a rich subject, ...

View answer
Non-conservative implementation implicit Euler
Accepted answer
13 votes

This might seem extreme, but this can be analysed exactly. Take the system $$ \dot x_1 = x_2, \qquad \dot x_2=-x_1, \qquad x_1(0) = 1, \qquad x_2(0)=0. $$ Let $X=(x_1,x_2)$ be the state vector, $\...

View answer
Complex Eigenvalues using eig (Matlab)
12 votes

You have an ill-conditioned eigenvalue problem. Consider a perturbation $\delta A$ in your matrix $A$—with double-precision floats this is $O(10^{-16})$. It turns the eigendecomposition from $$ X^{-1}...

View answer
Generate a set of orthogonal vectors to a given vector
Accepted answer
12 votes

There is a known mathematical question here: you are given a unit vector, lying on the $(n-1)$-sphere in $\mathbb{R}^n$, $v\in S^{n-1}$, and you would like to associate with each such vector a frame ...

View answer
Clenshaw-type recurrence for derivative of Chebyshev series
Accepted answer
11 votes

You can just take Clenshaw's recurrence $$ u_k(x) = 2xu_{k+1}(x)-u_{k+2}(x)+\color{red}{a_k},\\ f(x) = x u_1(x)-u_2(x)+\color{red}{a_0} $$ and differentiate it directly: $$ u_k'(x) = 2xu_{k+1}'(x)-u_{...

View answer
Is there any point to using hypot() for $\sqrt{1+c^2}$, $0 \le c \le 1$ for real numbers
Accepted answer
10 votes

Summary: it is cheapest and most accurate to use sqrt(fma(c, c, 1)) if you have FMA, and sqrt(1+c*c) otherwise. In my testing, though, the difference is extremely marginal: of the 1065353216 32-bit ...

View answer
Integrating Lagrange polynomials with many nodes, round-off
10 votes

You can evaluate this using the Björck-Pereyra algorithm for solving Vandermonde systems, because you are evaluating $b^\top V^{-1}$ with $b=(2,0,\frac23,0,\frac25,0,\ldots)$, and the algorithm is ...

View answer
Numerical evaluation of the Exponential Integral Ei by rational Chebyshev approximations fails
Accepted answer
10 votes

There is a prime missing both in your code and in the expression in your question. In the original paper, the expression is: $$ \log(x/x_0) + (x-x_0) \frac{\sum^\prime_{0\leq j\leq n}p_j T_j^*(x/6)}{\...

View answer
diagonalization of matrix - omitting small matrix elements
10 votes

It's not implementation-dependent in the sense that this is a mathematical operation performed on your matrix. However, it is very much matrix-dependent. If your matrix is diagonalizable and $A=XDX^{-...

View answer
Accurate computation of Gauss-Kuzmin entropy
Accepted answer
9 votes

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}{\...

View answer
Numerical calculation of Integral of Si(x)/x
Accepted answer
9 votes

For large $x$, you can use the same approach as the one in the paper you cite. They define these auxiliary functions that link $\mathrm{Si},\mathrm{Ci}$ with $\sin,\cos$, and the point of why this ...

View answer
Is it possible to proof a-b+b = a for all double floating-point numbers?
Accepted answer
9 votes

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

View answer
Nonlinear eigenvalue problem - MATLAB code
9 votes

Given a nonlinear eigenvalue problem of the form $A(\lambda)x = 0$, reducing it to a real equation $\det(A(\lambda))=0$ is known to be a poor method for just the reason you've discovered yourself. The ...

View answer
Methods for solving $x'=Ax+b$ for small, sparse, singular $A$
9 votes

Any general-purpose ODE solver should be able to handle this linear coupled system of ODE very easily, for example: scipy.integrate.ode CVODE from the Sundials solver suite; it appears to have Python ...

View answer
Name of an Optimization Approach to Reduce Size of Variable Space
Accepted answer
9 votes

This is coordinate descent. I believe it's used on very large-scale problems when other methods like gradient descent might be too slow (e.g., http://epubs.siam.org/doi/abs/10.1137/100802001). It ...

View answer
Why does LSODA fail to integrate the logistic function?
Accepted answer
8 votes

When you use $r=5$, the initial condition is $$ x(-10) \approx e^{-50} \approx 1.9\times 10^{-22}. $$ This is much smaller than the machine epsilon, $2\times 10^{-16}$, and it is very likely that ...

View answer
Stabilizing a 3x3 real symmetric matrix eigenvalue calculation
Accepted answer
8 votes

This is trying to compute the eigenvalues by computing the roots of the characteristic polynomial. In this case, the characteristic polynomial is $p(t) = t^3-2t^2x$, $x=1.25\times 10^6$, and zero is a ...

View answer
Numerically computing the advection equation
Accepted answer
8 votes

I see several issues: The DFT computed with fft puts the zero mode at the beginning of the array, and if you want to compute the derivative, it is necessary to apply fftshift/ifftshift to the array N ...

View answer
Finding the matrix inverse given a solver for the matrix equation $Ax=b$
Accepted answer
8 votes

Your two ideas make it much too complicated. If $X$ is the inverse of $A$, $$ AX=I, $$ and $x_i$ is the $i$-th column of $X$ and $e_i$ is the $i$-th column of the identity matrix $I$ ($e_i$ is a ...

View answer
Stability analysis of Heun's method
Accepted answer
8 votes

Notice that $\hat C_*=(1-r F(h))\hat C_n$, but the sign in front of $r$ is lost when you use $\hat C_*$ inside $\hat C_{n+1}$. I'm looking at p.149 of Numerical Solution of Time-Dependent Advection-...

View answer
Newton iteration for cube root without division
Accepted answer
8 votes

Cube roots are not nearly as important as square roots (e.g., for normalizing vectors), so that might be why they are discussed much less. In general, if you apply Newton's method to $x^\alpha-\beta$,...

View answer
The rate of convergence for finite difference methods for Poisson's equation with piecewise constant data
Accepted answer
8 votes

Let's examine the one-dimensional three-point stencil case in detail, because I think it's important to be clear just how this behaviour arises, and what it means to set a point to a certain value in ...

View answer
Numerical evaluation of an elliptic integral in python
Accepted answer
8 votes

The problem is almost definitely with how QUADPACK (which is the backend used by scipy.integrate.quad) handles numerically small integrands. Essentially the integrand is so small (at $x=0$ it is $6.58\...

View answer
C - OpenMP, MPI, Serial Program
Accepted answer
8 votes

I think some of your issues are more important than others and some of your emphasis is misplaced. In pursuing overhead, you are in danger of making your program unmaintainable. It is easier to write ...

View answer
Eigenvalues of $ab^T$
Accepted answer
7 votes

This is sometimes known as Buzano's inequality (http://www.jstor.org/stable/2159168). In general, if $\|x\|=1$, $P=xx^{\top}$ is a projection operator, so a simple application of Cauchy-Schwarz leads ...

View answer
1
2 3 4 5
7