83
votes
Accepted
What does "symplectic" mean in reference to numerical integrators, and does SciPy's odeint use them?
Let me start off with corrections. No, odeint doesn't have any symplectic integrators. No, symplectic integration doesn't mean conservation of energy.
What does ...
- 12.2k
44
votes
Accepted
Which Runge-Kutta method is more accurate: Dormand-Prince or Cash-Karp?
Since I just finished optimizing a lot of them in software, DifferentialEquations.jl, I decided to just lay out a comparison of the main Order 4/5 methods. The Fehlberg method was left out because it'...
- 12.2k
17
votes
What does "symplectic" mean in reference to numerical integrators, and does SciPy's odeint use them?
To complement Chris Rackauckas answer, to state some of the mathematical nonsense as well as some stuff you almost certainly know, a dynamical system is Hamiltonian if there is a description with ...
- 2,229
17
votes
Accepted
CUDA & Python for numerical integration and solving differential equations
Julia's DifferentialEquations.jl is all GPU-compatible. If you make your arrays GPU-based arrays, then the solver recompiles to be all on the GPU (no data transfers). For example:
...
- 12.2k
13
votes
Energy conservation in RK4 integration scheme in C++
RK4 is not symplectic so it has no guarantee of energy conservation. Especially when solving an N-body problem where two bodies pass by close to each other the energy conservation can be violated ...
- 2,391
13
votes
Accepted
How can I numerically integrate the Kepler problem?
I'm going to assume for the moment that your code is correctly implemented and that this problem isn't a bug.
Believe it or not, gradual increase of the energy is the expected behavior of most simple ...
- 10.1k
12
votes
Numerical evaluation of highly oscillatory integral
Use Plancherel's theorem to evaluate this integral.
The basic idea is that for two functions $f,g$,
$$
I=\int_{-\infty}^{\infty} f(x) g^*(x)dx = \int_{-\infty}^{\infty} F(k) G^*(k) dk
$$
where $F,G$...
- 663
10
votes
Integrating Lagrange polynomials with many nodes, round-off
You can evaluate this using the Björck-Pereyra algorithm for solving Vandermonde systems, because you are evaluating $b^\top V^{-1}$ with $b=(2,0,\frac23,0,\frac25,0,\ldots)$, and the algorithm is ...
- 11.4k
10
votes
Accepted
Why is velocity Verlet better than Verlet for gravity if it has a worse order of magnitude for the error term
That is a common misconception. Verlet in whatever form is a second-order method.
The misconception results from the fact that the truncation error of the Verlet formula is of 4th order. The naive ...
- 5,509
8
votes
Accepted
How to solve the problem without using symbolic computation
You can solve this numerically in Python without symbolic computation.
from __future__ import print_function, division
import numpy as np
from numpy import exp
from scipy.integrate import quad
from ...
7
votes
Accepted
Numerical integration of a hypergeometric function
First of all, from the first paragraph of your attempts at a solution, I assume that the $z_j$ are non-negative? In that case, the integrand has no real problematic points (it's monotonous, decreasing)...
- 6,169
7
votes
Accepted
What is wrong with this Euler method code in python?
I did not check your code, however, the result you are getting is also verified by scipy.integrate.odeint
...
- 731
7
votes
Accepted
How to estimate the error of trapezoidal rule using discrete data?
You may use the ideas of error extrapolation as one uses it to construct high-order Runge Kutta methods.
Depending on the function that you interpolate, the interpolation error $I - I_h$, where $I$ ...
- 3,418
7
votes
Accepted
Evaluating the surface integral in an FEM (Finite Elements Method) procedure
The particular surface integral you want to calculate is basically a specific case of integrating a function over a surface defined
in terms of two parametric coordinates. Lets first consider this ...
- 5,974
7
votes
Accepted
Integration of the Fermi distribution using Python
One of your problems is the system of units that you are using. Just changing the units improves the results
...
- 8,370
7
votes
Accepted
How do I integrate a function defined over an arbitrary area?
Instead of directly integrating over the area, it is often more convenient
to use the
divergence theorem to replace the area integral with an integral over the
boundary edges.
The divergence theorem ...
- 5,974
7
votes
Numerical evaluation of highly oscillatory integral
The key to the evaluation of oscillatory integrals is to truncate integral at the right point. For this example you need to choose upper limit of the form
$$
\pi\mathbb{N}+\frac{\pi}{2}
$$
Before ...
- 71
7
votes
Trapezoid rule vs Gaussian quadrature: what am I missing?
How accurate do you want the answer ?
How costly is evaluating your function ? If it is costly, then you dont want to use a rule with too many nodes.
How many times do you want to do the quadrature ?...
- 2,993
7
votes
Numerical integration giving incorrect sign
This might be an accuracy problem in computing the second term, because of those large exponentials when $x \gg 1$.
I would first work on that term: gather $e^x$ out from numerator and denominator and ...
- 11.1k
7
votes
Energy conservation in RK4 integration scheme in C++
Main error
As I pointed out in the previous questions of this series
RK4 integration of the three-bodies problem with C++
the primary problem is that the methods are not implemented correctly. The ...
- 5,509
6
votes
Accepted
Looking for an accurate algorithm to find the area of a oscillatory function
You can look up various quadratures. One method that should fair better is Gauss Quadrature.
I would also recommend looking into any adaptive quadrature schemes. There are many of them out there, so ...
- 3,898
6
votes
Integral in log-log space
You can just change variables. Setting $a=log(x)$, $b(a)=log(y(x))$. The integral becomes
$F(r)=\int^{log(r)}_{-\infty} exp(a+b) da$
You have to be a little careful because you are integrating from ...
- 483
6
votes
Radial integration of expensive function with Bessel weights
For the Hankel transform, one can classify the methods into four major groups:
Numerical quadrature-based.
Fourier-based ones.
Asymptotic expansion of Bessel into sines and cosines.
Projection-slice ...
- 8,602
6
votes
Accepted
How to solve this set of equations involving an integral?
It is always easier to solve a differential equation rather than an integral equation.
You can easily differentiate your last equation w.r.t the time variable $t$, and set the initial condition $\psi(...
- 1,628
6
votes
Accepted
Error on a integral quantity with noise
Assuming that you mean the following inequality in your prompt
$$ |F_\mathrm{true}(x) - F(x)| \le |dF(x)| \qquad \forall x,$$
a simple bound for $dM$ is the following
$$ dM \equiv | M_\mathrm{true} - ...
- 2,465
6
votes
Accepted
Evaluating an integral numerically at many points
The point of @WolfgangBangerth is exactly what I mentionend in my comment, so I'd always try this first. In the best case, with millions of partitions $[a_{i},a_{i+1}]_{i\in \{0,\ldots,N-1\}}$ (where $...
- 3,042
6
votes
Numerical integration in Python with unknown constant
What you want seems inherently impossible, and that’s not due to restrictions of Python.
The only way we can arrive at a situation where we only need to apply a single quadrature is to get ...
- 2,022
6
votes
Numerical evaluation of highly oscillatory integral
Ooura's method for Fourier sine integrals works here, see:
Ooura, Takuya, and Masatake Mori, A robust double exponential formula for Fourier-type integrals. Journal of computational and applied ...
- 2,125
5
votes
Accepted
Composite simpson's rule with odd intervals
A simple solution is to apply Simpson's (standard) rule to the first $n-3$ grid points, where $n-3$ is even for $n$ odd, and cover the remaining three gridpoints via the Simpson 3/8 formula:
$$I_{3/8}...
- 3,042
Only top scored, non community-wiki answers of a minimum length are eligible