We study the nonequilibrium time evolution of the Bose-Einstein condensate of interacting bosons confined in a leaky box, when its number fluctuation is initially suppressed. We take account of quantum fluctuations of all modes, including of the bosons. As the wave function of the ground state that has a definite number N of interacting bosons, we use a variational form which is obtained by operating a unitary operator on the number state of free bosons. Using we identify a “natural coordinate” of the interacting bosons, by which many physical properties can be simply described. The can be represented simply as a number state of we thus call it the “number state of interacting bosons” (NSIB). To simulate real systems, for which if one fixes N at N will fluctuate at later times because of a finite probability of exchanging bosons between the box and the environment, we evaluate the time evolution of the reduced density operator of the bosons in the box as a function of the leakage flux J. We concentrate on the most interesting and nontrivial time stage, i.e., the early time stage for which much earlier than the time when the system approaches the equilibrium state. It is shown that the time evolution can be described very simply as the evolution from a single NSIB at into a classical mixture, with a time-dependent distribution, of NSIBs of various values of N at Using we successfully define the cosine and sine operators for interacting many bosons, by which we can analyze the phase fluctuation in a fully quantum-mechanical manner. We define a new state called the “number-phase-squeezed state of interacting bosons” (NPIB), which is characterized by a complex parameter It is shown that for can be rewritten as the phase-randomized mixture (PRM) of NPIBs. Among many possible representations of this representation is particularly convenient for analyzing the phase fluctuations and the order parameter. We study the order parameter according to a few typical definitions, as well as their time evolution. It is shown that the off-diagonal long-range order (ODLRO) does not distinguish the NSIB and NPIB. Hence, the order parameter defined from ODLRO does not distinguish them, either. On the other hand, the other order parameter defined as the expectation value of the boson operator , has different values among these states. In particular, for each element of the PRM of NPIBs, we show that evolves from zero to a finite value very quickly. Namely, after the leakage of only two or three bosons, each element acquires a full, stable, and definite (nonfluctuating) value of