  We present a numerical method that consistently implements thermal
fluctuations and hydrodynamic interactions to the motion of Brownian particles
dispersed in incompressible host fluids. In this method, the thermal
fluctuations are introduced as random forces acting on the Brownian particles.
The hydrodynamic interactions are introduced by directly resolving the fluid
motions with the particle motion as a boundary condition to be satisfied. The
validity of the method has been examined carefully by comparing the present
numerical results with the fluctuation-dissipation theorem whose analytical
form is known for dispersions of a single spherical particle. Simulations are
then performed for more complicated systems, such as a dispersion composed of
many spherical particles and a single polymeric chain in a solvent.
