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

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

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

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

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

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

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

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

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