45

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 symplectic mean and when should you use it? First of all, what does symplectic mean? Symplectic means that the solution exists on a symplectic manifold. A symplectic manifold is a solution set ...


32

Since I just finished optimizing a lot of them in a 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's commonly known to be less efficient than the DP5 method. Backstories Dormand-Prince 4/5 The Dormand-Prince method was developed to be accurate as a 4/5 ...


14

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 coordinates $\mathbf{p}$ and $\mathbf{q}$ and a functional, $\mathcal{H(\mathbf{p},\mathbf{q})}$ such that $$\frac{d\mathbf{q}}{dt}=+\frac{\partial \mathcal{H}}{\...


11

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$ are the Fourier transforms of $f,g$. Your functions both have relatively small support in the spectral domain. Here, $\sin x / x \rightarrow \text{rect}(k)$ ...


9

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 known to be forward-stable (see Error analysis of the Björck-Pereyra algorithms for solving Vandermonde systems by Nick Higham, http://www.maths.manchester.ac.uk/~...


8

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 scipy.optimize import root def f1(a1, a2, x): return exp(a1 * x + a2 * x * x * x) / (1 + x * x) def f2(a1, a2, x): return exp(a1 * x + a2 * x * x * x) *...


8

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$ is the actual integral value and $I_h$ the value obtained by the piecewise trapezoidal rule, may have a smooth asymptotic expansion $$ I - I_h = C_p h^p + C_{p+...


7

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 say that it is mostly confined within $[a-\sigma,a+\sigma]$, then split the integral as $$ \int_l^u f(x) \; dx = \int_l^{a-\sigma} f(x) \; dx + \int_{a-\...


7

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 general case. A fundamental relation from differential geometry is (using your notation) $$ d{\bf \Gamma_e} = \left\{\begin{array}{c} \frac{\partial x}{\partial\...


7

One of your problems is the system of units that you are using. Just changing the units improves the results import numpy as np import scipy.integrate as integrate eigenvalue = [0.9, 1.3] fermi = 1.0 T = 300 kB = 8.6173303e-5 def fermi_integral(E, fermi, T): return 1 / (1 + np.exp((E - fermi) / (kB * T))) for i in range(len(eigenvalue)): result = ...


7

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 explaining why it should work, let me at first show that it actually produces good results. Asymptotics It is easy to guess that asymptotic series has the ...


6

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 searching up algorithms for adaptive integration would help you. This wikipedia link should give you some insight into adaptive integration.


6

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 $-\infty$. What you have to do exactly will depend on what $y(x)$ looks like.


6

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 behavior you should expect. If you use a more accurate integration scheme you will see faster convergence. You can easily implement the trapezoid rule for ...


6

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). The only difficulty is the infinite integration range. And then it all depends on the accuracy you want... Do you want machine epsilon precision (or there ...


6

I did not check your code, however, the result you are getting is also verified by scipy.integrate.odeint from scipy.integrate import odeint def ode(x, t, N, theta, mu, k): dxdt = (N - x) * mu - theta * x / (k + x**2) return dxdt N = 50 theta= 0.6 mu = 0.6/412.5 k = 1 x0 = 0 t = np.linspace(0, 30, 101) sol = odeint(ode, x0, t, args=(N, ...


6

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(0) = \psi_0$ for the following differential equation: $$\frac{d\psi}{dt}=\psi(1+\kappa \ln{(\psi)}+x_0\psi).$$ Once you have the solution for the above, ...


6

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} - M | = \left|\int_{x_1}^{x_2} F_\mathrm{true}(x) - F(x) dx\right| \le \int_{x_1}^{x_2} |dF(x)| dx. $$ This also assumes that $dF$ is absolute integrable.


6

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 $a_{i} < a_{i+1}$), the length of these is quite small such that some low-order method could be appropriate. This will require roughly $N \times (s-1)$ ...


6

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 in three dimensions for a vector function $F$ can be written $$ \int_V \nabla \cdot {\bf F} dv = \int_S {\bf F} \cdot ds $$ That is, an integral over the ...


6

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 analytically get rid of all dependencies of $T$ in the integral. To this end, the best we can do is to apply the substitution $x=Ty$ to your integral: $$ I(T) = \int_0^∞ \...


6

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 mathematics 112.1-2 (1999): 229-241. I wrote an implementation of this algorithm but never put in the work to get it fast (by, say caching nodes/weights), but ...


5

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 play with the time step of an integrator (e.g., making it smaller by a factor of ten) and you get a different result, then the time step was too large; if the ...


5

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 seems extremely likely that the error is on your end when you're modifying the time step. If it is possible (i.e. not too computationally expensive). I would ...


5

For 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 methods. The following paper gives a nice overview of this methods (types of methods): M. J. Cree and P. J. Bones, "Algorithms to numerically evaluate the Hankel ...


5

Note: I'm somewhat worried at this point that the integral values Mathematica gives me are bogus. I thought it was working because it gave a sensible-looking result in a short time, but it might be the case that the method it tries to use is buggy or that I did something wrong. So it might be that the code below isn't working at all, I don't know, sorry. ...


5

I want to mention another idea in case it helps. The truncation error of replacing $\int_{-\infty}^{\infty}$ with $\int_{-L}^L$ is on the order of $$\begin{aligned} \int_L^\infty \cos(tx)g(x)\,\mathrm{d}x &= \frac{1}{t}\sin(tx)g(x)\big|_{L}^{\infty} - \int_L^\infty\frac{\sin tx}{t}g'(x)\,\mathrm{d}x \\&= -\frac1L g(L)\sin(tL) + \text{asymptotically ...


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

What are some general heuristics for using "lsoda" vs. "dopri5"? I suppose an obvious one consideration: "is your problem stiff?". Well, if the problem I am working with is stiff, then "lsoda" is struggling with "Excess work done." too...so, the BDF method doesn't seem to help too much with stiffness. If forced to choose between only those two ...


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