Acoustic wave propagation in a one-dimensional waveguide connected with Helmholtz resonators is studied numerically. Finite amplitude waves and viscous boundary layers are considered. The model consists of two coupled evolution equations: a nonlinear PDE describing nonlinear acoustic waves, and a linear ODE describing the oscillations in the Helmholtz resonators. The thermal and viscous losses in the tube and in the necks of the resonators are modeled by fractional derivatives. A diffusive representation is followed: the convolution kernels are replaced by a finite number of memory variables that satisfy local ordinary differential equations. A splitting method is then applied to the evolution equations: their propagative part is solved using a standard TVD scheme for hyperbolic equations, whereas their diffusive part is solved exactly. Various strategies are examined to compute the coefficients of the diffusive representation; finally, an optimization method is preferred to the usual quadrature rules. The numerical model is validated by comparisons with exact solutions. The properties of the full nonlinear solutions are investigated numerically. In particular, the existence of acoustic solitary waves is confirmed.