  We have developed an efficacious algorithm for simulation of the beam-beam
interaction in synchrotron colliders based on the nonlinear $\delta f$ method,
where $\delta f$ is the much smaller deviation of the beam distribution from
the slowly evolving main distribution $f_0$. In the presence of damping and
quantum fluctuations of synchrotron radiation it has been shown that the slowly
evolving part of the distribution function satisfies a Fokker-Planck equation.
Its solution has been obtained in terms of a beam envelope function and an
amplitude of the distribution, which satisfy a coupled system of ordinary
differential equations. A numerical algorithm suited for direct code
implementation of the evolving distributions for both $\delta f$ and $f_0$ has
been developed. Explicit expressions for the dynamical weights of
macro-particles for $\delta f$ as well as an expression for the slowly changing
$f_0$ have been obtained.
