We develop and implement higher order methods for the derivative free solution of shape optimizationproblems. As model problem we consider a generalized version of Bernoulli's exterior free boundary problem which involves a source term and non-constant boundary data. We derive first and second order shape update rules based on the asymptotic expansion of the solution of either the pure Dirichlet boundary value problem or the mixed boundary value problem. In particular, besides the pointwise boundary update, also variational boundary updates will be considered.Analogies to the Hadamard shape gradient and shape Hessian will be studied in order to prove convergencerates of the optimization scheme.

The solution of the state equation is computed by a Nystr\"om-based $p$-boundary element method using trigonometric polynomials. It enables in combination with a parametric representation of the free boundary the computation of higher order derivatives of the solutionat the boundary. The developed shape optimization algorithm will be extended to solve constrained shape optimization problems like the electromagnetic shaping or the optimization of elastic bars.