People who code: we want your input. Take the Survey
13

A fairly simple method would be to choose a basis in function space and convert the integral transformation to a matrix. Then you can just invert the matrix. Mathematically, here's how that works: you need some set of orthonormal basis functions $T_i(x)$. (You can get away without them being normalized too, but it's easier to explain this way.) Orthonormal ...


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


8

The usual approach to testing optimization algorithms is to compare how many function evaluations they need to find the minimum (or get within a fixed tolerance $\varepsilon$ of the minimum). This is easily implemented in any code, and it does not matter in that case how fast your computer is, or what the resolution of the tic/toc mechanism is. There are a ...


6

Arnold Neumaier's answer is discussed in more detail in section 3.2 of the paper "Approximating Spectral Densities of Large Matrices" by Lin Lin, Yousef Saad and Chao Yang (2016). Some other methods are also discussed but the numerical analysis at the end of the paper shows that the Lanczos method outperforms these alternatives.


6

I hope I understood the question correctly. They try to compute exactly the same thing, so they really are equivalent. I'll use Chebyshev polynomials because they are easy to analyze. Given a function $f(x)$ on $[-1,1]$, the spectral interpolant is the truncation of $$ \begin{aligned} f(x) &= \sum_{n\geq0} \bar a_n T_n(x), \\ \bar a_n &= \frac{1+[n&...


5

It may look like the error is getting larger, but I think you are confusing multiple issues here. The first is that you are looking at finite approximating Taylor polynomials $p_3(x), p_5(x), p_7(x)$ evaluated at $x=\pi/2$. Taylor's theorem says that for analytic functions (which the sine is), $p_N(\pi/2)\rightarrow 1$ as $N\to \infty$. Your results appear ...


5

You could try using a library that implements the Fast Multipole Method (FMM), which should drastically reduce the amount of memory you need and will decrease the complexity of matrix-vector products from $\mathcal{O}(N^{2})$ to $\mathcal{O}(N)$. It is difficult to implement, but there should be some libraries out there. Another algorithm for N-body ...


5

The "fast" version is (damped/over-relaxed) Jacobi instead of SOR, which is a multiplicative method. Even with damping factor of 1.0, Jacobi is not guaranteed to converge, as you would see if you applied it to a 1D problem, or added extreme anisotropy to your multi-dimensional problem. Your SOR ("slow") implementation is slow because it is written in pure ...


5

In order to inconvenience as many people as possible, long ago, mathematicians and physicists decided to use two different conventions on whether $\theta$ or $\phi$ is the latitude angle. Greengard's notes use the physicists' convention that $\theta$ is latitude and $\phi$ is longitude, whereas Scipy uses $\theta$ for longitude and $\phi$ for latitude, so ...


5

You can write a curve as a parametric equation ($x(t)$, $y(t)$) for a range say $0 < t < 1$. Then you can have a Fourier decomposition of both $x(t)$ and $y(t)$. A rectangle can be written as such a parametric equation. However, you want to decompose a given function as a sum of functions from a specific class. There is a whole field devoted to this. ...


5

I believe that for odd polynomial degrees you can get superconvergence by one order at the jump. For even degrees I don't think that's the case. For the odd case, we consider a smooth function $u:\mathbb{R}\to\mathbb{R}$. Consider two adjacent intervals $I^-=[x_{i-1},x_i]$ and $I^+=[x_i,x_{i+1}]$ both of length $h$, and denote the $L^2$ projection of $u$ on ...


5

I think this is a really cool question, which might have a really cool answer, if someone was willing to think about it hard enough to publish on. But when we're talking smooth functions which need to be evaluated quickly, I reach for B-splines. Create the pairs $\{x_j, f(x_j)\}_{j=1}^{n}$, and use (say) a cubic B-spline (here's an implementation for the ...


4

Taylor series is not a good way to do this; it takes a lot of terms to get reasonable accuracy. I answered a similar question on stackoverflow a while ago. Here's the body of the answer: Here are some good slides on how to do power series approximations (NOT Taylor series though) of trig functions: http://www.research.scea.com/gdc2003/fast-math-...


4

Randomized algorithms can accurately approximate your matrix if it has rapidly decaying singular values (i.e. a Gaussian sampling of the range is likely to pick up the action of the matrix). Since your matrix is symmetric, you could approximate a truncated eigenvalue decomposition using random linear algebra and simply take the square root of the eigenvalues....


4

Let's use the example (approximate the square root function on $[0.25,1.0]$ with a quartic polynomial) to step through your calculations. I suspect that the code is going to work with only modest changes. Having chosen an initial set of six interpolation points: $$ x_k = 0.25, 0.4, 0.55, 0.7, 0.85, 1.0\;\; (k=0,\ldots,5) $$ we proceed to interpolate the ...


4

Fast filter approximations have been studied for a long time, especially to implement IIR filters, like Gaussians and their derivatives. You may want to reuse such concepts, with keywords like integral image, summed-area tables, box filter, recursive filtering. You can start from a recent review: A Survey of Gaussian Convolution Algorithms, 2013 and dig ...


3

Your confusion mainly seems to be with calculus, not with numerical methods. Specifically, make sure you understand the difference between an indefinite integral and a definite integral. consider what the initial conditions or boundary conditions are for your problem. Perhaps you know $q(t_0) = q_0$ for a given time $t_0$. Write the integral formulation ...


3

I will answer for the simplest case: $\frac{dy}{dt} = -\lambda{y}, \hspace{4mm} y(0)=1$ Note: You need an initial condition which you did not specify in your original question. In this simple problem $q(t) = y(t)$ and $p(t) = -\lambda{q(t)}$. If you use a basic forward Euler finite difference we get: $\frac{y^{n+1}-y^{n}}{dt} = -\lambda{y^{n}}, \hspace{...


3

It seems you are looking for an optimization (here, specifically a fit). You have some experimental data and you want to have a simple model explaining it. The simplest way is using the Ordinary Least Squares. The idea is to minimize the standard deviation or the error you make on the vertical axis. Another one would be the PCA where you try to find an ...


3

There are two standard choices for the initial approximation of $B$ in BFGS: Either you choose $B=\frac{\|g_0\|}{\delta} I$ where $\|g_0\|$ is the gradient in the very first iteration and $\delta$ a "typical step size" from $x_k$ to $x_{k+1}$. Or you choose $B=\frac{y_1^T y_1}{y_1^T s_1}I$ using the standard notation used for BFGS. The (excellent) book by ...


3

In the 1-dimensional case, there are three main methods which are similar to what you are describing which come to my mind. The secant method uses the derivative approximation $$f'(x_n)\approx\frac{f(x_n)-f(x_{n-1})}{x_n-x_{n-1}}$$ to compute the next iteration. As André Nicolas points out, the secant method is quite desirable for its fast convergence and ...


3

Disclaimer: I'm new(er) to 'R' and so I don't know the sparse FFT library (SE,MIT,Berk) off the top of my head. This stuff has been published for ~2 years. I would be surprised if it doesn't exist - it should and would be valuable addition to the language if it doesn't. If I were highly motivated (not at the moment), I would make it myself. The question ...


3

It's a bit heavy, but not all of it must be read to get a useable method : Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow, Bao and Du. I believe -- it's been a while since I worked on it -- that you have to move to imaginary time (as you said, by using $-it$ instead of $t$). Then you want to propagate along ...


3

By using orthonormal polynomial basis functions, then you do not need to solve a linear set of equations in order to get the polynomial coefficients, similar to the Fourier series coefficients. Because for large $N$ the matrix, of the linear set of equations, will get ill-conditioned. It is at least required that the basis functions are orthogonal, but ...


3

This approach is the result of a Taylor approximation, which assumes that we are not very far from the solution (solution is the case when $R\approx I$ - no more updates can be made). Under large rotations this linearization will fail. Moreover, point to plane distance is prone to sliding errors when normals and points have particular configurations, such ...


3

The question you ask is simpler than you may think. For any mesh size $h$, you can compute $\lambda(h)$. If this dependence is continuous, then you can do a Taylor expansion of $\lambda(h)$ around $h=0$ to obtain exactly the form you want. The only question in your case is whether $\lambda(h)$ is a continuous function. This is made slightly more ...


3

As with most questions about the computation of special functions, the Digital Librarary of Mathematical Functions is a good place to start. In particular, see chapter 6, which deals with the exponential integral and $\mbox{li}(x)$. You'll find that different methods (e.g. the power series for small $x$ versus asymptotic expansion for large $x$) work best ...


3

In a one-step method $$ y_{n+1}=y_n+h\Phi_f(x_n,y_n,h) $$ one gets a truncation error for the exact solution $$ y(x_{n+1})=y(x_n)+h\Phi_f(x_n,y(x_n),h)+h^{p+1}\tau(x_n) $$ For the error propagation of $e_n=y_n-y(x_n)$ this gives approximately $$ e_{n+1}=e_n+h\partial_y\Phi_f(x_n,y(x_n),h)e_n-h^{p+1}\tau(x_n)+... $$ The error is then of order $p$, $e_n=c(x_n)...


2

I found a paper (circa 2005) by Boyar et al which studies algorithms for problems of this type, which they call maximum resource bin-packing problems (to contrast with the usual goal of bin-packing to minimize a resource such as the number of bins). Usually, for bin packing problems, we try to minimize the number of bins used or in the case of the dual ...


2

You can reduce the problem to computing the singular value decomposition, for which there exist many fast methods and codes. For fast 3x3 SVD, I found this paper. To reduce the polar decomposition to the SVD, suppose the polar decomposition is written in the following form, $$M = U P,$$ with orthogonal $U$ and positive semidefinite $P$ (Ie., $Q \...


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