Polynomial approximations of curves This is the 3D version of this question. The responses to that question contained a lot of complaints about fuzzy definition of the problem, so I made this new question very narrow and explicit.
For a given fixed degree, $n$, I'm interested in constructing parametric polynomial approximations of degree $n$ of certain curves, for example those that arise as intersections of two surfaces. This is important in engineering and manufacturing computing, where parametric polynomials are usually the only curves available (in the form of Bézier curves).
A typical example is as follows:
Let $f(x,y,z) = x^2 + y^2 - 16$, and let $S(f)$ be the cylindrical surface where $f(x,y,z)=0$. Similarly, let $g(x,y,z) = y^2 + z^2 - 25$, and let $S(g)$ be the cylindrical surface where $g(x,y,z)=0$. Let's focus on the curve $C$ that is the portion of the intersection of these two surfaces lying in the positive octant where $x,y,z \ge 0$. A few simple calculations show that $C$ has end-points $\mathbf{p}_0 = (4,0,5)$ and $\mathbf{p}_1 = (0,4,3)$.

I want to construct a parametric polynomial curve $t \mapsto \mathbf{x}(t) = \big(x(t), y(t), z(t)\big)$ of degree $n$ that approximates $C$ and such that $\mathbf{x}(0) = \mathbf{p}_0$ and $\mathbf{x}(1) = \mathbf{p}_1$. The error in the approximation will be measured by
$$
E(\mathbf{x}) = \sup\big\{  f(\mathbf{x}(t))^2 + g(\mathbf{x}(t))^2 : 0 \le t \le 1\big\}
$$
First question: for a given $n$ how can I find the polynomial curve $\mathbf{x}$ of degree $n$ that minimizes $E(\mathbf{x})$? 
If we were approximating a real-valued function, there would be a great deal of approximation machinery that would help us: the Weierstrass approximation theorem, the equi-oscillation criterion, the Remez algorithm, etc. If I parameterize $C$, then each of its three components will be a real-valued function, and I can apply these known techniques. But the result will depend on how I parameterize $C$, and will not be optimal. 
Second question: Has any of the above-mentioned approximation machinery been generalized beyond the real-valued case?
 A: Do you care more about the curve or the surfaces? Replacing the first function in your example by $100(x^2 + y^2 - 16)$ or $x^2-z^2+9$ gives the same curve but a different error function. 
Also, consider the intersection of the cylinder $(x-11)^2+(y-11)^2=221$ and the plane $x+y-z=1.$ It is an ellipse and the portion $\mathit{C}$ in the positive octant goes from $(1,0,0)$ to $(0,1,0)$ visiting $(k,k,2k-1)$ along the way where $k=11+\sqrt{442}/2 \approx 21.51.$ I'm guessing that a polynomial curve minimizing what you want would go along the small missing section which goes through $(j,j,2j-1)$ where $j= 11-\sqrt{442}/2 \approx 0.49$ so $2j-1<0$. Of course we could specify that the entire curve should stay in the positive octant but then $(t,1-t,0)$ might be pretty hard to beat with anything we would think of as approximating $\mathit{C}.$. 
That said, your particular example seems well behaved. In the case of $n=2$ (chosen for illustration) you seek $9$ constants so that for $\mathbf{x}(t)=(at^2+bt+c,dt^2+et+f,gt^2+ht+i)$ the maximum of $ f(\mathbf{x}(t))^2 + g(\mathbf{x}(t))^2$ on $[0,1]$ is minimized. The extra condition on the position for $t=0$ and $t=1$ reduces to $3$ degrees of freedom. 
In the case of polynomial approximation of a real curve $f(t)$ by a degree $n$ polynomial $p_n(t)$ on $[0,1]$ (with or without side conditions for $t=0,1$) It is usually not possible to explicitly find the  best polynomial approximation although it is known that there is one which is unique and it can be recognized by the extreme values of $f-p_n$ in $[0,1]$ being equal (say to $E$) in absolute value, alternating in sign, and sufficient in number. (In particular, $f-p_n=0$ (at least) $n+1$ times in $[0,1].$) There are iterative schemes , as I recall,  which can do well. Pick  a set of points $X=\{t_1,t_2,\cdots,t_n\}$ in the open interval (somehow) and solve the system $p_n(0)=f(0),p_n(1)=f(1)$ and $p_n(t_i)-f(t_i)=(-1)^iE$ for the unknown coefficients of $p_n$ and the unknown error value $E.$ Find the set $Y$ of maximum errors of $|p_n-f|.$ If by a miracle you got the right thing you can recognize it. Otherwise, replace $X$ with $Y$ and try again. I'm not sure how to adapt that.
I was able to get a pretty good approximation for your example with $n=2$ $$(- 3.776\,{t}^{2}+ 7.776\,t,- 2.831\,{t}^{2}- 1.169\,t+4,- 0.813\,{t}^
{2}+ 2.813\,t+3).
$$ The worst errors are about $1.1$ near $t=0.2$ and $t=0.8.$ This arose as follows:
Specify that the error is $0$ (as required) at $t=0,1$ and also at $t=1/2$ This gives one degree of freedom (Say for $d$ which corresponds to  picking which point on the curve corresponds to $t=1/2$.) The error increases  to a maximum in the first half then decreases to $0$ at $t=1/2$ then increases then decreases again. Keep adjusting until the two maxima are roughly equal.
LATER Thinking more about your particular example, you have a curve $\mathit{C}$ so that the projections onto the $xy$ and $zy$ planes are (quarter) circles. For a curve with  projections onto two orthogonal planes which are easily described curves your method with the intersection of the two corresponding (generalized) right cylindrical surfaces seems like a natural choice based on the curve itself. Even better might be projection (if possible) onto three pairwise orthogonal planes. In your case the projection onto the $xz$ plane is is a portion of the hyperbola $z^2-x^2=9$. So minimize the maximum value over the interval of $ f(\mathbf{x}(t))^2 + g(\mathbf{x}(t))^2+h(\mathbf{x}(t))^2.$ Of course another three orthogonal planes would give a different error function though maybe not too different a result. Still, for an application, Bézier curves for each component might be fast and satisfactory.
