The governing equations of a two-phase flow have a strong nonlinear term due to the interactions between gas and water such as capillary pressure, water saturation, and gas solubility. This nonlinearity is usually ignored or approximated in order to obtain analytical solutions. The impact of such ignorance on the accuracy of solutions has not been clear so far. This study seeks analytical solutions without ignoring this nonlinear term. Firstly, a nonlinear mathematical model is developed for the two-phase flow of gas and water during shale gas production. This model also considers the effects of gas solubility in water. Then, iterative analytical solutions for pore pressures and production rates of gas and water are derived by the combination of travelling wave and variational iteration methods. Thirdly, the convergence and accuracy of the solutions are checked through history matching of two sets of gas production data: a China shale gas reservoir and a horizontal Barnett shale well. Finally, the effects of the nonlinear term, shale gas solubility, and entry capillary pressure on the shale gas production rate are investigated. It is found that these iterative analytical solutions can be convergent within 2-3 iterations. The solutions can well describe the production rates of both gas and water. The nonlinear term can significantly affect the forecast of shale gas production in both the short term and the long term. Entry capillary pressure and shale gas solubility in water can also affect shale gas production rates of shale gas and water. These analytical solutions can be used for the fast calculation of the production rates of both shale gas and water in the two-phase flow stage.