# Tag Info

65

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

35

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

17

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: using OrdinaryDiffEq, CUDA, LinearAlgebra u0 = cu(rand(1000)) A = cu(randn(1000,1000)) f(du,u,p,t) = mul!(du,A,u) prob = ODEProblem(f,u0,(0.0f0,1.0f0)) # Float32 is better on ...

15

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}}{\... 12 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) ... 10 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 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 ... 7 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, ... 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

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

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

7

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 ? If it has to be repeated many times, like in finite element methods, the cost of computing weights/nodes of gauss rule will not matter much since it has to be ...

7

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 simplify it out (leaving a summand $e^{-x}$ in the numerator). However, rather than guessing it might be necessary to understand better where the accuracy ...

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

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

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

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 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 the Clenshaw-Curtis quadrature weights w^{\text{cc}}_k on the Chebyshev extrema grid y_k for k=0,\ldots,n:$$ y_k = \cos\left(\frac{k\pi}{n}\right)\\ w^{\...

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