The Mumford-Shah model has been well acknowledged as an important method for image segmentation. This paper discussed the problem of simultaneous image segmentation and smoothing by approaching the Mumford-Shah paradigm from a numerical approximation perspective. In particular, a novel iterative relaxation algorithm for the numerical solving of the Mumford-Shah model was proposed. First, the paper presented mathematically the existence of a solution in the weak formulation of GSBV space. Second, some approximations and numerical methods for computing the weak solution were discussed. Finally, a minimization method based on a quasi-Newton algorithm was put forward. The proposed algorithm found accurately the absolute minimum of the functional at each iteration. Considering the important role of a discrete finite element approximation method in the sense of Γ-convergence, an adjustment scheme for adaptive triangulation was applied to improve the efficiency of iteration. Experimental results on noisy synthetic and jacquard images demonstrate the efficacy of the proposed algorithm.