I derived the following equation in my research but I don't understand the equation very well. I hope to get some help here.

$$\displaystyle \frac{\partial f(r,t)}{\partial t} = h_0(t) f(r,t) \Big( \bar{r}(t) - r \Big)$$, where

• $$r \ge 0$$ is the physical quantity of interest
• $$\displaystyle f(r,t)$$ is the probability density function of $$r$$ at time $$t$$
• $$h_0(t)$$ is a given continuous nonnegative function: $$h_0(t) \ge 0$$ for $$\forall t \in \mathbb{R}$$
• $$\displaystyle \bar{r}(t) \equiv \int_0^{\infty} r \, f(r,t) dr$$  is the mean of $$r$$ at time t

My questions:

1. Has anyone seen any equations like this one? What kind of equation is it? And can we say anything about its properties, in addition to the ones listed below?
2. Any comments on the strategies to analytically solve the equation, given some initial condition $$f(r, 0)$$ (eg, log-normal distribution with parameters $$\mu$$ and $$\sigma$$)?
3. To numerically solve the equation given $$f(r, 0)$$, currently I represent $$f(r,t)$$ as a vector over a fixed grid of $$r$$, and update it with Euler steps $$f(r,t+dt) = f(r,t) + h_0(t) f(r,t) \big( \bar{r}(t) - r \big) dt$$. It works, but seems a bit slow, and unstable: for some large $$dt$$ and $$r \gg \bar{r}(t)$$, $$f(r, t+dt)$$ can become negative, as opposed to staying very close to zero. Any comments on how to improve the numerical integration algorithm here?
4. I am ultimately only interested in $$\bar{r}(t)$$, and am solving $$f(r,t)$$ just for $$\bar{r}(t)$$. Given the information presented here, is there a way to solve for $$\bar{r}(t)$$ directly?

Many thanks!

So far I have discovered the following properties of the equation:

1. Conservation of probability mass: $$\displaystyle \int_0^{\infty} f(r, t) dr = 1 \implies \int_0^{\infty} f(r, t+dt) dr = 1$$
2. Decay of mean: At time $$t$$, for $$r < \bar{r}(t)$$ its density increases because $$\displaystyle\frac{\partial f(r,t)}{\partial t} \propto (\bar{r}(t) - r)>0$$, while the opposite is true for $$r > \bar{r}(t)$$. As a result, the distribution of $$r$$ shifts to the left over time, towards a Dirac's delta function at $$r=0$$. In fact, one can derive that $$\displaystyle \frac{d \bar{r}(t)}{dt} = - h_0(t) V(t)$$, where $$\displaystyle V(t) \equiv \int_0^{\infty} \big( r - \bar{r}(t) \big)^2 f(r,t) dr$$ is the variance of $$r$$ at time $$t$$. That is, the mean of $$r$$ decays at a rate proportional to its variance.
3. Invariance of Gamma distribution under the evolutionary operator: If $$f(r,0)$$ is Gamma distributed with shape parameter $$k$$ and scale parameter $$\theta_0$$, then $$f(r,t)$$ is also Gamma distributed with parameter $$k$$ and $$\displaystyle \theta(t)=\frac{\theta_0}{\theta_0 \int_0^t h_0(\tau) d\tau +1}$$.