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


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

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

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

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

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

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


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


5

Indeed the problem was that you were trying to calculate your nonlinear term incorrectly and you forgot that: $$\mathcal{F}[(f(x))^{2}] \neq (\mathcal{F}[f(x)])^{2}$$ I changed your code to this: import numpy as np import matplotlib.pyplot as plt from matplotlib.pyplot import cm nu = 1 L = 100 nx = 1024 t0 = 0 tN = 200 dt = 0.05 nt = int((tN - t0) / 0....


5

The discrete Fourier transform for a signal of period $T$ with $N$ samples reads in its inverse or reconstruction form as $$ y(t)=\frac1{N}\sum_{k=-N/2}^{N/2}c_k e^{i2\pi k\frac{t}{T}} $$ with redundancy in $c_N$ and $c_{-N}$ if $N$ is even. Sampling this at points $t_m=\frac{mT}{N}$ gives a completely determined linear system whose solution is given by the ...


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

It is so trivial to pick a solution on a box domain for the Laplace equation. Just pick a function, say $\bar u(x,y,z)=x^2y^2\sin(z)$ (chosen in a way so that it isn't in your ansatz space), then compute $f=-\Delta \bar u$ and solve the Laplace equation $$ -\Delta u = f. $$ Of course, your numerical solution should converge to $\bar u$. If you can't deal ...


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

When I change your code to this: psi[x_] := Piecewise[ {{0, x < -2} , {CubicLagrangeInterpolation[{-3, -2, -1, 0}, {0, 0, 0, 1}, x], -2 <= x <= -1} , {CubicLagrangeInterpolation[{-2, -1, 0, 1}, {0, 0, 1, 0}, x], -1 <= x <= 0} , {CubicLagrangeInterpolation[{-1, 0, 1, 2}, {0, 1, 0, 0}, x], 0 <= x <= 1} , {...


3

It's a question of function spaces. If your function is in $L^2$, then you know that you can approximate the function using a Fourier series because the Fourier basis (i.e., the sines and cosines) are dense in $L^2$. In other words, you know that if you take more and more terms of the series, then you will approximate the function better and better. On the ...


3

I haven't done the algebra, but I think you have the formula wrong. http://mathworld.wolfram.com/FourierTransformGaussian.html. There's some factors of $\pi$ missing from your answer. Which scaling of the Fourier Transform did you use?


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


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


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