  Simulation of the dynamics of dust-gas circumstellar discs is crucial in
understanding the mechanisms of planet formation. The dynamics of small grains
in the disc is stiffly coupled to the gas, while the dynamics of grown solids
is decoupled. Moreover, in some parts of the disc the concentration of the dust
is low (dust to gas mass ratio is about 0.01), while in other parts it can be
much higher. These factors place high requirements on the numerical methods for
disc simulations. In particular, when gas and dust are simulated with two
different fluids, explicit methods require very small timestep (must be less
than dust stopping time $t_{\rm stop}$ during which the velocity of a solid
particle is equalized with respect to the gas velocity) to obtain solution,
while some implicit methods requires high temporal resolution to obtain
acceptable accuracy. Moreover, recent studies underlined that for Smoothed
particle hydrodynamics (SPH) when the gas and the dust are simulated with
different sets of particles only high spatial resolution $h<c_{\rm s} t_{\rm
stop}$ guaranties suppression of numerical overdissipation due to gas and dust
interaction.
  To address these problems, we developed a fast algorithm based on the ideas
of (1) implicit integration of linear (Epstein) drag and (2) exact conservation
of local linear momentum. We derived formulas for monodisperse dust-gas in
two-fluid SPH and tested the new method on problems with known analytical
solutions. We found that our method is a promising alternative for the
previously developed two-fluid SPH scheme in case of stiff linear drag thanks
to the fact that spatial resolution condition $h<c_{\rm s} t_{\rm stop}$ is not
required anymore for accurate results.
