22

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


15

Computing the term $e^x$ is definitely significantly more expensive than computing a lower-order polynomial -- say $x^4$. But it may be ten to 100 times more expensive at most, not "crazy" expensive. So I suspect that if Matlab takes forever to compute something, then that's because the character of your ODE changes significantly. For example, the presence ...


13

The singular value decomposition for a symmetric matrix $A=A^{T}$ is one and the same as its canonical eigendecomposition (i.e. with an orthonormal matrix-of-eigenvectors), while the same thing for a nonsymmetric matrix $M=U \Sigma V^T$ is just the canonical eigenvalue decomposition for the symmetric matrix $$ H=\begin{bmatrix}0 & M\\ M^{T} & 0 \end{...


11

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


11

The problem with binomial coefficients is that they are a quotient of factorials, which may be rather large and in particular much larger than the final binomial coefficent: $$\binom{n}{k} = \frac{n!}{k!(n-k)!}$$ For example, $\binom{10}{5}=252$ but $10!=3628800$ (the numerator). In another example $\binom{100}{5}=75287520$ and thus can be represented as a ...


10

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)}{\sum^\prime_{0\leq j\leq n}q_j T_j^*(x/6)}. $$ This primed sum is very common and standard when dealing with Chebyshev series: it means that the first term of ...


9

This is more of a question appropriate to the CS forums, but here is the short of it: it requires you to know how the compiler parses numbers. In particular, when the compiler sees something like -2147483648, it parses this as unary minus applied to 2147483648. But the latter (a positive number) is too large for the int data type, and so may be wrapped ...


8

A possible solution (suggested by @gammatester) is to use Jacobi polynomials. This circumvents the problem of catastrophic cancellation in adding the large polynomial coefficients by 'naive' polynomial evaluation. The radial Zernike polynomial can be expressed by Jacobi polynomials as follows (see equation (6)) $$ R^m_n(\rho) = (-1)^{(n-m)/2}\rho^m \cdot P^...


7

Intermediate results may overflow even if the final result doesn't. There are ways one can re-arrange an equation so that this doesn't happen but it's up to the programmer (maybe an optimizing compiler that's also a theorem solver can do it too - but I've never seen one) and the CPU won't re-order a sequence of arithmetic instructions on its own. Instead of ...


7

In this paper, Honarvar and Paramesran derive an interesting method to compute the radial Zernike polynomials in a very nice recursive way. The recursion formula is surprisingly straightforward, without division or multiplication by large integers: $$ R^m_n(\rho) = \rho \left(R^{|m-1|}_{n-1}(\rho)+R^{m+1}_{n-1}(\rho)\right) - R^{m}_{n-2}(\rho)$$ I'd ...


5

I think there's a simple way to do this. You have a rational function of identical cosh/sinh terms, where every expression is a homogeneous polynomial in cosh/sinh, and the only problem is that these exponential terms overflow. The function does not diverge as these terms approach infinity, so if you divide every numerator and denominator by the same power ...


5

No it's not OK. There are no negative literals per-se in C. "-12345" is interpreted as the combination of the positive literal 12345 with the unary minus operator. 0x80000000 is outside the range of int and will be interpreted as a literal of type unsigned int. On most compilers you could get around this by using a typecast #define INT_MIN ((int)(...


4

While $a^{b}$ is usually computed via $\exp(b \log (a))$, in real-life high-quality implementations of functions like pow() the logarithm is computed with extra precision to counter the error magnification effect of $\exp$, as demonstrated in this answer for example. This suggests that we should favor computation via a built-in general exponentiation ...


4

It's not overcomplicated: That's how $x^y$ is actually defined -- via $$ x^y = \exp(y \log x). $$ The point is that for non-integer $y$, it's not at all obvious what $x^y$ actually means: It is not just a product of $x$ with itself, repeated $y$ times. The definition of $x^y$ then is done by pulling the expression back to functions we have previously ...


4

If we want to approximate $$f(x) := 1 - \frac{\sqrt{1+x^{2}}}{1+\frac{1}{2}x^2}$$ on $[10^{-10}, 1]$ as accurately as possible, and in particular with an accuracy close to the machine precision associated with a given IEEE-754 floating-point format, neither naive computation nor a Taylor series expansion seem suitable. There is an issue with subtractive ...


4

I took an example from your other question and tried to demonstrate here what is going on. So, for completeness: We are applying the 6th order finite-difference scheme for numerical differentiation using the following code. First, let's apply it on the constant function of different scale ($C_1, C_2, C_3$): $$ U(x)=C,\quad x\in[0,1] $$ Since $\frac{\text{d}...


2

This a followup to my comment about smoothing the expression for $k$ that is a function of $x$ and $\dot x$. I converted your expression for $k$ into the following function to make experimentation easier. The input variable xxdot is just $x\dot x$ function k=calcK(xxdot) k = 4; if(xxdot>0) k = k + 4; end end Then I wrote a second version of ...


2

Kirill's suggestion to use symbolic linear algebra is on-point. Another approach is to work out the characteristic polynomial over a number of finite fields, then use the Chinese remainder theorem to reconstruct the characteristic polynomial over $\mathbb{Z}$. Finite field linear algebra does not suffer from the same kinds of numerical problems as you'll ...


2

If, as you say, you only really need the characteristic polynomial of your matrices, and not the eigenvalues (solving for the roots of the characteristic polynomial doesn't in general give accurate eigenvalues), and since your matrices are really quite small ($40\times40$) and have only integers, it might be sufficient to just use matlab's symbolic ...


2

If you do not want to (or just can't) juggle with the expression, just use the multiprecision library of your choice. Below an example with Boost.Multiprecision in C++. This has the invaluable advantage that you do not have to touch the original expression at all, thanks to operator and function overloading. #include <iostream> #include <cmath>...


1

To tackle your issue of very, very small numbers, you need to use a arbitrary precision arithmetic library, like MPFR. https://www.mpfr.org/ MPFR is awesome and will continually jack up the precision until it is sufficient to avoid roundoff error or you run out of memory. In my experience, I never used more than a 128 bit mantissa (a double has something ...


1

You can easily approximate $f(x)$ up to fourth order accuracy by using Taylor expansion as: $$f(x) = 1 - \frac{\sqrt{1+x^{2}}}{1+\frac{x^{2}}{2}}$$ Taylor expansions: $$\sqrt{1+x^{2}} = 1 + \frac{x^{2}}{2} - \frac{x^{4}}{8} + O(x^{5})$$ $$\frac{1}{1+\frac{x^{2}}{2}} = 1 - \frac{x^{2}}{2} + \frac{x^{4}}{4} + O(x^{6})$$ So: $$f(x) = 1 - \frac{\sqrt{1+x^{...


1

Although the question has a great answer, here's a rule of thumb for small singular values, with a plot. If a singular value is nonzero but very small, you should define its reciprocal to be zero, since its apparent value is probably an artifact of roundoff error, not a meaningful number. A plausible answer to the question "how small is small ?" ...


1

I'm trying to answer it with physical intuition and arguments. I assume you'll get $p(\mathbf{r},t)$ as solution. Take only one dimension for example i.e. writing as $p(x,t)=\Sigma_n \Sigma_m a_{mn} Sin(k_mx-\omega_nt)$. If you Fourier decompose it in space first, you'll $k_m$s. Then for each $k_m$, do Fourier transform in time to get $\omega_n$. Now plot ...


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