Abstract The problem of approximating a given set of data points by splines composed of Pythagorean hodograph (PH) curves is addressed. We discuss this problem in a framework that is not only restricted to PH spline curves, but can be applied to more general representations of shapes. In order to solve the highly non-linear curve fitting problem, we formulate an evolution process within the family of PH spline curves. This process generates a family of curves which depends on a time-like variable t. The best approximant is shown to be a stationary point of this evolution process, which is described by a differential equation. Solving it numerically by Euler's method is shown to be related to Gauss–Newton iterations. Different ways of constructing suitable initial positions for the evolution are suggested.