Skip to main content
85 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 ...
Chris Rackauckas's user avatar
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'...
Chris Rackauckas's user avatar
18 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 ...
origimbo's user avatar
  • 2,249
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: ...
Chris Rackauckas's user avatar
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 ...
helloworld922's user avatar
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 ...
Daniel Shapero's user avatar
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$...
smh's user avatar
  • 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 ...
Kirill's user avatar
  • 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 ...
Lutz Lehmann's user avatar
  • 6,064
10 votes

Implementation of Monte-Carlo Integration

From Wikipedia: The naive Monte Carlo approach is to sample points uniformly on Ω[...] There is an implicit assumption here that a uniform distribution on $\Omega$ exists. It is well-known that such ...
whpowell96's user avatar
  • 2,443
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 ...
numerical-solution's user avatar
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 ...
Stelios's user avatar
  • 731
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)...
GertVdE's user avatar
  • 6,179
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$ ...
Jan's user avatar
  • 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 ...
Bill Greene's user avatar
  • 6,064
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 ...
nicoguaro's user avatar
  • 8,510
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 ...
Bill Greene's user avatar
  • 6,064
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 ...
David Saykin's user avatar
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 ?...
cfdlab's user avatar
  • 3,028
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 ...
Federico Poloni's user avatar
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 ...
Lutz Lehmann's user avatar
  • 6,064
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 ...
spektr's user avatar
  • 4,248
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 ...
Anton Menshov's user avatar
  • 8,672
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(...
HBR's user avatar
  • 1,638
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} - ...
Richard Zhang's user avatar
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 $...
davidhigh's user avatar
  • 3,127
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 ...
Wrzlprmft's user avatar
  • 2,032
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 ...
user14717's user avatar
  • 2,155
5 votes
Accepted

Integrating Lagrange polynomials with many nodes, round-off

The calculation of $$ \int_{-1}^{1} L_k(x)\,\text{d} x $$ for the Lagrange polynomials $L_k$ defined on an arbitrary grid $x_k, k=0,\ldots,n$ can be performed by the following two steps: Calculate ...
davidhigh's user avatar
  • 3,127

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