15

For integrations with many variables, the Monte Carlo method usually is a decent fit. Its error decreases as $O (\sqrt{N})$ where N is the number of equidistributed points selected. Of course this is not good for low dimension (1D and 2D) spaces where high order methods exist. Most of these deterministic methods, however, take a large number of points in ...


12

Your matrix $A$ isn't a circulant matrix- it's just Toeplitz. The method that you're trying to use fundamentally only works for circulant systems. Furthermore, your $a$ vector doesn't have the "-1" in it anywhere, so you clearly don't have sufficient information. A method that involves embedding the $n$ by $n$ Toeplitz matrix in a double-sized ...


11

Have you thought of using Barycentric Interpolation? A detailed description on how to do it efficiently for Chebyshev nodes is given in Section 5 of this paper. This is actually an exact evaluation of the Chebyshev interpolant. If you're evaluating a polynomial of degree $n$ at $m$ nodes, the cost is in $\mathcal O(nm)$. Update Another alternative, if you ...


10

Let me first answer all the questions: What is the theoretical convergence rate for an FFT Poison solver? The theoretical convergence is exponential as long as the solution is sufficiently smooth. How fast should this energy converge? The Hartree energy $E_H$ should converge exponentially for a sufficiently smooth solution. If the solution is less ...


9

17. Lots and lots of work has gone into good fft implementations and it's unlikely you'll be able to reliably outperform a good fft library. For example, fftw "automatically adapts itself to your machine, your cache, the size of your memory, the number of registers, and all the other factors that normally make it impossible to optimize a program for more ...


9

This is a linear integral equation (because the right hand side is linear in $h$) and presumably the problem is ill-posed (for one, because of the multiplication by $\sin(\theta)$ but also because the "integral kernel" is smooth). Methods to solve this (approximately) are general Galerkin methods or collocation methods. For Galerkin methods you take a ...


8

Sparse grid quadrature is an alternative approach to integrate in higher dimensions. Quadrature relies on evaluating a weighted sum of function values at specific "optimal" points. Traditional quadrature uses a tensor product grid construction in higher dimensions, which means that you would have to evaluate the function at an exponentially growing ...


8

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 to make sure the derivative is correct. It is easy to see for yourself what the correct expression is by working it out with pen and paper, and see also the ...


7

Michael Pippig at the University of Chemnitz (Germany) has implemented an MPI-parallelized FFT that uses FFTW in the background. This might help you: http://www-user.tu-chemnitz.de/~mpip/software.php?lang=en It is using the algorithm proposed by Plimpton from Sandia National Labs as suggested by Eldila's comment.


7

The frequency content of the interpolated signal is significantly influenced by the interpolation basis. If you have a band-limited function that you have adequately sampled (i.e. satisfying Nyquist criterion), interpolating with any function that is not band-limited to the same frequency will indeed introduce high-frequency noise. Unfortunately, exact band-...


6

In a word, yes. A least squares problem arises in the context of an over determined system. That means that you have more equations than unknowns and would be the case if you had more data points than Fourier coefficients that you're trying to find. This would give you the "best fit" Fourier series of a given order. For example, you could find the best fit ...


6

This turned out much longer than I planned so I hope it is useful and can be a resource for others with the same questions. I believe you mean to ask primarily about the discrete Fourier transform (in which case your definition (2) for $\hat{f}_\alpha$ is not right, for the DFT it should be $\hat{f}_\alpha = \sum_{j=0}^{2M-1} x_j \exp(- i\alpha\frac{\pi j}{...


6

You are solving Poisson's equation: $$\nabla^2 \Phi = 4\pi G \rho.$$ Notice that if $\Phi$ is a solution, then so is $\Phi+C$ for any constant $C$. Furthermore, $C$ will have no effect on your simulation since the forces depend only on derivatives of $\Phi$. The divide-by-zero issue is just a manifestation of this degeneracy. In the Fourier-...


6

From your link, consider the definition of the round off error and the statement "Since the exact solution must satisfy the discretized equation exactly, the error must also satisfy the discretized equation." This is actually only true if the PDE is homogeneous, that is, if we can write it in the form $\mathcal{L}(u;u_t,u_x,u_{xx},\ldots)=0$, with all terms ...


6

Running your code, it seems like your pulse looks kinda like this: (sorry for not adding units to the plots, I used the same as you, i.e. t is in fs and w in rad/fs) So, the FWHM is not correct (should be 4 fs) and the angular carrier/central frequency is messed up (should be ca. 2.3 / ca. 0.37 per fs). I think there is two things not quite right here: 1) ...


5

This is pure haar scaling function approximation of $f(t)$ $$f(t)=\sum_{k=-\infty}^{\infty}a_k\phi(2^jt-k)$$ where $$\phi(t) = \begin{cases}1 \quad; 0 \leq t < 1,\\0\quad;\mbox{otherwise.}\end{cases}$$ is haar scaling function, $j$ scaling aparameter, $k$ translation parameter and the coefficients $$a_k=\int f(t)\phi(2^jt-k)dt$$ So on our case if $f(t)...


5

Yes, it is true, the CVX package (disclosure: I am the author) does not support the use of the fft() command in constraint and objective expressions. However, it's relatively simple to get the equivalent result, in both the 1D and 2D cases. First, 1D. Given a vector $x\in\mathbb{R}^n$, the DFT of $x$ is equal to $Wx$, where $W$ is the so-called DFT matrix. ...


5

The easiest way to generate a divergence-free 2D function is to use a streamfunction, $\psi$. Where $$ u=\frac{\partial\psi}{\partial y}\quad \rm{and}\quad v=-\frac{\partial \psi}{\partial x}$$ Once you generate a random $\psi$, you can then use your Fourier representation (or finite differences) to compute the two derivatives. I haven't tried this kind of ...


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

1) You have misinterpreted the order the fft2 output elements are stored. This is a very common mistake, since it is very easy to get confused. Thankfully, there is an easy way to avoid such bugs. Arrange the Fourier space vectors the "natural way" % Fourier space vectors kx = 2*pi/Lx*[-Nx/2:Nx/2-1]'; ky = 2*pi/Ly*[-Ny/2:Ny/2-1]'; Then, make sure the fft2 ...


5

The point is to guarantee that the Fourier coefficients exist, not to make really fine distinctions about what might happen if the guarantee is violated. Furthermore, you would be hard pressed to find a reasonable function for which this is an issue, if it exists at all. The class of integrable functions is already so wide as to include everything you might ...


4

As an alternative you could look into using the Goertzel Algorithm to directly compute the frequency components you are interested in.


4

The FFT can be used for periodic boundary conditions. Because von Neumann boundary conditions are effectively "mirror" boundary conditions, you have to do a "mirrored continuation", before you can apply an FFT. One drawback of this approach is that you will blow up the data volume by a factor 4 (which is not important if you are just interested in ...


4

I want to expand on my comment and rework the example you reference in a way that should be more understandable than the original and to explain why fft returns the coefficients the way it does. For reference, the fft portion of the example is : Nx = size(x,2); k = 2*pi/(b-a)*[0:Nx/2-1 0 -Nx/2+1:-1]; dFdx = ifft(1i*k.*fft(f)); d2Fdx2 = ifft(-k.^2.*fft(f)); ...


4

The interpolation indeed affects the Fourier transform. @Steve already gives the correct answer in general, but I want to give you an example that helps the intuition more. Think for example that you have samples a sine function on a set of equally spaced points. If you did a (discrete) FFT on these points alone, you would of course recover a Fourier ...


4

I did some more investigating about the accuracy of the standard recurrence relation method vs. the newer Bunck algorithm. It seems that in fact the Bunck algorithm is generally more accurate for all $n$ and $x$. I made a python script using the mpmath arbitrary precision library to calculate the Hermite functions to lots of decimal places, and then compared ...


3

Ah, I have realized the answer to my own question: It is important to recognize that the initial data $v_0,...,v_N$ is not stored on a uniform grid, but rather at the Chebyshev points $$ x_j = \cos\frac{\pi j}{N},\qquad j=0,...,N.$$ Now as long as the initial data has a decent polynomial interpolation, then \begin{align} v_j = p(x_j) &= a_0 + a_1x_j + \...


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