Differentiating the solutions to differential equations

A kind of problem that turns up quite often in physics and elsewhere is to find the solution of an ordinary differential equation, given some initial conditions. That is, you have a system with some state $x$ represented as a vector of real numbers, you have an initial state $x_0$, and you have a rule describing the evolution of the state: $$ \frac{dx}{dt} = F(x,t) $$ And your goal is to find $x(t)$. This standard problem has some standard solution techniques, some quite advanced - Runge-Kutta methods, symplectic methods, Hermite integrators. A few are implemented in scipy. But it sometimes happens that solving this problem is only part of a more complicated process, say of fitting, where it would be nice to have the derivatives of the solution with respect to the various initial conditions. It turns out this isn't too hard to work out, usually.

Full post

Weighted Poisson Uncertainties

I recently ran across a rather awkward mathematical problem. I'm trying to make a histogram of photon arrival phases, complete with an uncertainty on the number of photons in each bin. Normally this is done by just taking the square root of the number of photons, which is at least approximately right based on Poisson statistics. But in my problem — data from the Fermi space telescope, which is not very good at localizing low-energy gamma rays — the photons are weighted: for each photon I have a probability that it really came from the source. So the values in the histogram should be the total probability. But what should the uncertainty be? The short version is: the square root of the sum of the squares of the weights.

ETA: This is in the literature, without justification as far as I can tell. See below.

Full post