19

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


13

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


9

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


8

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


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


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

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