Abstract.A new numerical method based on fictitious domain methods for shape optimization problems governed by the Poisson equation is proposed. The basic idea is to combine the boundary variation technique, in which the mesh is moving during the optimization, and efficient fictitious domain preconditioning in the solution of the (adjoint) state equations. Neumann boundary value problems are solved using an algebraic fictitious domain method. A mixed formulation based on boundary Lagrange multipliers is used for Dirichlet boundary problems and the resulting saddle-point problems are preconditioned with block diagonal fictitious domain preconditioners. Under given assumptions on the meshes, these preconditioners are shown to be optimal with respect to the condition number. The numerical experiments demonstrate the efficiency of the proposed approaches.Mathematics Subject Classification. 49M29, 65F10, 65K10, 65N30, 65N55.