76 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 ...
user avatar
39 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'...
user avatar
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 ...
user avatar
  • 2,199
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: ...
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$...
user avatar
  • 653
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 ...
user avatar
  • 11.4k
8 votes

Numerical integration of sharp peaked function (position of peak known)?

If you know where the peak is, then you can always split the interval. For example, if you know that the peak is at $a$ and has a "width" (however you want to define that) of $\sigma$ so that you can ...
user avatar
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 ...
user avatar
8 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$ ...
user avatar
  • 3,408
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)...
user avatar
  • 6,091
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 ...
user avatar
  • 731
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 ...
user avatar
  • 5,744
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 ...
user avatar
  • 8,111
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 ...
user avatar
  • 5,744
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 ...
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 ?...
user avatar
  • 2,983
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 ...
user avatar
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 ...
user avatar
  • 3,683
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 ...
user avatar
6 votes
Accepted

Simple integration in Matlab is off

You're using a first order accurate integration technique (the rectangle rule) and your error is proportional to $1/N\propto\Delta r,\Delta \phi,\Delta \theta$. This is exactly the type of convergence ...
user avatar
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 ...
user avatar
  • 8,452
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(...
user avatar
  • 1,606
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} - ...
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 $...
user avatar
  • 2,792
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 ...
user avatar
  • 1,914
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 ...
user avatar
  • 2,021
5 votes

Numerical integration using RKF7(8) - different results

If you get qualitatively different results from two ODE integrators, then the time step choice of at least one of them is too large. Which one that is is not immediately obvious to say, but if you ...
user avatar
5 votes

Numerical integration using RKF7(8) - different results

I think the first step is to confirm which of the solutions is more accurate. If you are using a reference implementation of RKF7(8) that has presumably been validated by others on other problems it ...
user avatar
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 ...
user avatar
  • 2,792

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