Japan J. Indust. Appl. Math., 13 (1996), 495-517 An Approximate Factorization/Least Squares Solution Method for a Mixed Finite Element Approximation of the Cahn-Hilliard Equation E.J. DEAN, R. GLOWINSKI and D.A. TREVAS Department of Mathematics, University of Houston, Houston, Texas 77204-3476, USA Received June 26, 1995 We discuss in this axticle the numerical solution of the Cahn-Hilliard equation mod- elling the spinodal decomposition of binary alloys. The numerical methodology combines a second-order finite difference time discretization with a mixed finite element space ap- proximation and a least squares formulation based on an approximate factorization of a fottrth-order elliptic operator which appears in the numerical model. The least squares problem - which is linear - is so]ved by a preconditioned conjugate gradient algorithm. The results of numerical experiments illustrate the possibilities of the methods discussed in this article. Key words: Cahn-Hilliard equation, spinodal decomposition, approximate factorization, mixed finite element approximations, preconditioned conjugate gradient algorithms 1. Introduction Immiscible binary mixtures such as Fe-A1 alloys enjoy under certain circum- stances (cooling below a critical temperature, for example) a phase separation phe- nomenon known as spinodal decomposition. Starting from the pioneering work of J.W. Cahn and J.E. Hilliard (refs. [1], [2]) this phase transition phenomenon has motivated a large number of investigations, from various points of view: modelling (refs. [1]-[4]), mathematical (refs. [5], [6], [7, pp. 147-158]) and computational (refs. [5], [8]-[16]). It is generally agreed that spinodal decompositions of binary mixtures are modelled by the Cahn-HiUiard equations, completed by appropriate initial and boundary conditions. Suppose that ~2 C ~I~ d (d = 1, 2, or 3) is the space region in which the spinodal decomposition is taking place; a prototype Cahn-Hilliard equation is given by: Ou (1.1) Ot A(U3--U)-}-s in ~2X (0, T) (1.2) u(x, 0) = ~0(x), = e ~, Ou 0 (1.3) On onAU = 0 on 0s In (1.1)-(1.3), (0, T) is the time interval on which the spinodal decomposition d is considered (with 0 < T _< +oo); 0/2 is the boundary of ~2; x = {x~}~=1 is a generic