# Tag Info

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

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

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

7

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

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

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

Regarding point 1: generally, applying a discretized derivative (stencil) twice is not equivalent to applying an equivalent (i.e. same number of points) stencil for the second derivative. Example in matrix representation: A three point stencil for the first derivative would correspond to a tridiagonal matrix. The same is true for the three point stencil for ...

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

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

1) One usually chooses one or two representative wall-normal distances and presents the spectra for <u'u'>, <v'v'>, etc. over the homogeneous directions. For example, see Figure 11 within Spalart 1988 (http://dx.doi.org/10.1017/s0022112088000345). 2) By "scale" do you mean a lengthscale? If so, a lengthscale quantifying what aspect of the flow?...

3

I asked this question on the Signal Processing exchange as well. Eventually I found the answer myself and posted it there. Here I copy/paste my answer as posted there. After many weeks I give the answer to my own question. There is a limit in which we can solve this problem in a reasonably simple way. Suppose we sum enough points in our DFT that the ...

3

1) If you want code, use Matthew Emmett's PyWENO package. Otherwise, the reference you give is the answer. You just need to do some simple calculus; you won't find that written down in a paper. 2) That kind of analysis does not apply directly to nonlinear schemes like WENO. If you can come up with a meaningful generalization of it, you should publish ...

3

Addition to the accepted answer. Here is a link to an open source implementation of Greengard's and Lee's method: https://finufft.readthedocs.io/en/latest/ It has wrappers to C, fortran, MATLAB, octave, and python. I believe the FINUFFT is written in C++. It is maintained and used at NYU Courant institute, SFU, Flatiron institute (obviously), University of ...

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