A cubic Bézier curve can replace a smooth function on a finite interval with four points, two of which are fixed by the function values. The approximation is compact, fast to evaluate, and compatible with vector graphics and geometry kernels. Finding a convincing curve, however, requires more than placing two controls by eye. The endpoint tangents must be preserved, the parametric curve must remain a graph, and the word best must refer to a precisely stated error measure.
Select a function below. The red line is the function, the blue line is its fitted cubic, and the dashed line is the control polygon. Every example is fitted in the browser from the same two-parameter model; none uses stored control points.
The Cubic Model
Let a continuously differentiable function be given on \([a,b]\). Its cubic Bézier approximation is the parametric curve
\[ \mathbf B(t)=(1-t)^3\mathbf P_0+3t(1-t)^2\mathbf P_1 +3t^2(1-t)\mathbf P_2+t^3\mathbf P_3, \qquad 0\leq t\leq1. \]
The factors \3\ in the two inner Bernstein terms are essential. Omitting them destroys the partition of unity, changes the endpoint values, and no longer describes the standard cubic Bézier curve. Writing the controls as two-dimensional points makes another important fact explicit: \t\ is a curve parameter, not the function's \x\-coordinate.
Interpolation at the interval endpoints fixes two controls immediately,
\[ \mathbf P_0=(a,f(a)),\qquad \mathbf P_3=(b,f(b)). \]
A direct count would now leave four scalar coordinates in \\mathbf P_1\ and \\mathbf P_2\. Endpoint derivatives reduce these to two much more meaningful unknowns.
Endpoint Tangents as Directions
Differentiating the cubic at its endpoints gives
\[ \mathbf B'(0)=3(\mathbf P_1-\mathbf P_0),\qquad \mathbf B'(1)=3(\mathbf P_3-\mathbf P_2). \]
The graph of \f\ has tangent direction \(1,f'(x))\ wherever the slope is finite. Normalize these directions,
\[ \mathbf d_a=\frac{(1,f'(a))}{\sqrt{1+f'(a)^2}},\qquad \mathbf d_b=\frac{(1,f'(b))}{\sqrt{1+f'(b)^2}}, \]
and place the internal controls along their endpoint tangent rays:
\[ \boxed{ \mathbf P_1=\mathbf P_0+\alpha\mathbf d_a, \qquad \mathbf P_2=\mathbf P_3-\beta\mathbf d_b, } \qquad \alpha,\beta\geq0. \]
The original eight control coordinates have become two nonnegative handle lengths. More importantly, this formulation uses directions rather than slope equations. A vertical tangent needs no division and no exchange of the coordinate axes: use \\mathbf d=(0,1)\ or \(0,-1)\. The square-root example in the figure starts with exactly such a tangent.
The scale of a direction matters. Unit vectors make \\alpha\ and \\beta\ geometric distances, which is useful for bounds and diagnostics. Using the unnormalized vectors \(1,f')\ is mathematically equivalent, but changes the numerical values assigned to the two parameters.
What Does Best Mean?
A Bézier curve is parametric: \\mathbf B(t)=(x(t),y(t))\. Comparing \y(t)\ with \f(a+(b-a)t)\ silently assumes a linear relationship between \t\ and \x\, even though moving the controls also changes \x(t)\. The vertical graph residual must instead be evaluated at the point's actual horizontal coordinate:
\[ r(t;\alpha,\beta)=y(t;\alpha,\beta)-f\bigl(x(t;\alpha,\beta)\bigr). \]
Several legitimate objectives follow from this residual. Uniform weighting in parameter space gives
\[ E_t(\alpha,\beta)=\int_0^1 r(t;\alpha,\beta)^2\,dt. \]
For a function graph, it is usually more natural to give equal horizontal intervals equal influence. If \x(t)\ is monotone, the substitution \dx=x'(t)\,dt\ yields
\[ \boxed{ E_x(\alpha,\beta)=\frac{1}{b-a} \int_0^1 r(t;\alpha,\beta)^2x'(t;\alpha,\beta)\,dt. } \]
The interactive fitter minimizes \E_x\ by composite Simpson quadrature. Its reported root-mean-square error is \\operatorname{RMSE}=\sqrt{E_x}\. MSE and RMS have the same minimizer; a mean absolute error or a maximum error generally does not. The maximum vertical error
\[ E_\infty=\max_{0\leq t\leq1}|r(t)| \]
is often the better acceptance criterion for rendering, but its nonsmooth minimax problem requires a different solver. Orthogonal distance, radial error, and curvature error define still other approximations.
A Metric Changes the Optimum
Consider the quarter circle \f(x)=\sqrt{1-x^2}\ on \[0,1]\. The familiar symmetric cubic uses equal handle lengths near \0.55228475\ because it is designed around circular geometry. Minimizing the horizontal-uniform vertical error above instead gives approximately
\[ \alpha=0.55418,\qquad \beta=0.54972. \]
The unequal result is not a defect: vertical distance distinguishes the horizontal and vertical axes. A cubic circle approximation optimized for radial or geometric deviation answers a different question. No set of controls is universally best without naming the metric.
Remaining a Function Graph
The residual \f(x(t))\ is meaningful even when the curve loops, but a looping curve is not a single-valued approximation of \f\. The required condition is
\[ x'(t)>0\quad\text{for }0<t<1. \]
Because \x'(t)\ is a quadratic Bézier polynomial, its minimum can be checked exactly from its endpoints and interior stationary point. A simpler conservative constraint is to require ordered control abscissae,
\[ P_{0x}\leq P_{1x}\leq P_{2x}\leq P_{3x}. \]
This sufficient condition is what the interactive implementation enforces. It permits a zero horizontal derivative at a vertical endpoint while preventing a reversal in the interior. Rejecting infeasible controls is preferable to taking an absolute value of \x'(t)\, which would hide the loop and optimize the wrong object.
Solving the Two-Parameter Problem
Substituting the tangent controls into the Bézier equation makes \x(t)\ and \y(t)\ affine in \\alpha\ and \\beta\. The residual is nevertheless nonlinear because \f\ is evaluated at the moving coordinate \x(t)\. Its derivatives are especially simple. With
\[ b_1(t)=3t(1-t)^2,\qquad b_2(t)=3t^2(1-t), \]
the chain rule gives
\[ \frac{\partial r}{\partial\alpha} =b_1(t)\left(d_{a,y}-f'(x(t))d_{a,x}\right), \]
\[ \frac{\partial r}{\partial\beta} =-b_2(t)\left(d_{b,y}-f'(x(t))d_{b,x}\right). \]
Sampling these residuals and derivatives produces a two-column Jacobian. The Gauss-Newton method, preferably with Levenberg-Marquardt damping, is therefore a natural solver for large batches of related fits. Bounds on the handle lengths and a line search are still needed to preserve monotonicity.
For the small interactive problem, a bounded Nelder-Mead simplex is a pragmatic choice. It needs no numerical derivatives, works directly with the quadrature objective and its infeasible boundary, and is inexpensive in only two dimensions. Several deterministic starting scales reduce dependence on a single initial simplex. This does not prove that every nonlinear objective has a unique global minimum; the returned controls are the best feasible solution found from those starts.
function fitFunctionGraph(f, derivative, a, b) {
P0 = [a, f(a)]
P3 = [b, f(b)]
da = unitTangent(derivative(a))
db = unitTangent(derivative(b))
objective(alpha, beta) {
P1 = P0 + alpha * da
P2 = P3 - beta * db
reject unless P0.x <= P1.x <= P2.x <= P3.x
return simpsonIntegral(residual(t)^2 * dxdt(t)) / (b - a)
}
return boundedMultiStartMinimize(objective)
} When One Cubic Is Not Enough
Four endpoint conditions leave only two shape parameters. That compression is the point, but it also limits what one segment can reproduce. A single cubic is a poor model when the interval contains several inflections, sharp changes of curvature, nearly singular derivatives, or oscillations. A low RMS can also conceal a narrow local excursion, so both RMS and maximum error should be checked.
Adaptive piecewise fitting gives a dependable remedy:
- Fit one cubic on the current interval.
- Evaluate the maximum graph error using more points than the optimization quadrature.
- Accept the segment when the tolerance is met.
- Otherwise split near the largest error or a known critical point and fit both subintervals.
- Reuse the true derivative at each split so adjacent segments share a tangent direction.
Matching tangent directions gives geometric \G^1\ continuity. Matching the complete derivative vectors requires coordinating the neighboring handle lengths and parameter scales, producing parametric \C^1\ continuity. If curvature continuity matters, the two segments must be optimized together with additional constraints.
Practical Checks
- Verify that function values are finite at both endpoints.
- Represent vertical endpoint tangents as vectors, not enormous substitute slopes.
- State the error norm and weighting whenever control points are published.
- Check \x'(t)\ over the complete interval, not only at sampled parameters.
- Use multiple initial guesses or a bounded global search when local minima are plausible.
- Measure the final curve with denser, independent samples and a maximum-error check.
- Split at singularities, discontinuities, and intervals whose error tolerance cannot be met.
The essential modeling decision is now visible. Endpoint values and tangent directions determine the admissible family; two handle lengths select a member of that family; the chosen metric decides which member deserves to be called best. Keeping those three layers separate turns an informal drawing technique into a reproducible approximation method.
References
- [Farin2002]Gerald Farin, Curves and Surfaces for CAGD: A Practical Guide, fifth edition, Morgan Kaufmann, 2002.
- [Nocedal2006]Jorge Nocedal and Stephen J. Wright, Numerical Optimization, second edition, Springer, 2006.
- [Nelder1965]John A. Nelder and Roger Mead, “A Simplex Method for Function Minimization,” The Computer Journal, volume 7, issue 4, 1965, pages 308–313.