Cubic Bezier curve tangent to Circle

Question

I have a cubic Bezier curve (4 points: P0, P1, P2, P3) : $$B(t)=((1-t)^3P_0+3t(1-t)^2P_1+3t^2(1-t)P_2+t^3P_3) , t \in [0,1]$$

I have a circle ("half-bottom-circle") defined as: $(a, b)$ its center, $R$ its radius, and parametric form : $$\left(u+a\ ,\ -\sqrt{R^{2}-u^{2}}+b\right), u \in [-R,R]$$

I would like to define $y_2$ (y-coordinate of $P_2$), so that when I move $y_1$, the circle and the Bezier are ALWAYS tangent. For them to be tangent, 2 conditions are necessary:

1- The slopes of the tangents (derivatives) at the point of tangency are equal : $$\frac{3\left(1-t\right)^{2}\left(y_{1}-y_{0}\right)+6\left(1-t\right)t\left(y_{2}-y_{1}\right)+3t^{2}\left(y_{3}-y_{2}\right)}{3\left(1-t\right)^{2}\left(x_{1}-x_{0}\right)+6\left(1-t\right)t\left(x_{2}-x_{1}\right)+3t^{2}\left(x_{3}-x_{2}\right)}=\frac{u}{\sqrt{R^{2}-u^{2}}}$$

2- The "y" coordinates at the point of tangency are equal : $$\left(1-t\right)^{3}y_{0}+3\left(1-t\right)^{2}ty_{1}+3\left(1-t\right)t^{2}y_{2}+t^{3}y_{3}=-\sqrt{R^{2}-u^{2}}+b$$

Now I have two unknowns "$t$" and "$u$", and I am looking for $y_2$. Algebraically I think it's going to be complicated. I am willing to use Newton-Raphson, but I don't know how to move forward. A little help would be welcome. Thank's. See my graph

Graph

Answer 1

Given the cubic Bézier curve of parametric equations:

$$ \begin{aligned} & b_x(u) := x_0\,(1 - u)^3 + 3\,x_1\,u\,(1 - u)^2 + 3\,x_2\,u^2\,(1 - u) + x_3\,u^3 \\ & b_y(u) := y_0\,(1 - u)^3 + 3\,y_1\,u\,(1 - u)^2 + 3\,y_2\,u^2\,(1 - u) + y_3\,u^3 \\ \end{aligned} \quad \quad \text{with} \; u \in [0,\,1] $$

and the circle of parametric equations:

$$ \begin{aligned} & c_x(v) := x_C + r\,\cos v \\ & c_y(v) := y_C + r\,\sin v \\ \end{aligned} \quad \quad \text{with} \; v \in [0,\,2\pi) $$

in order to be tangent it must be:

$$ \begin{cases} b_x(u) = c_x(v) \\ b_y(u) = c_y(v) \\ b_x'(u) = \lambda\,c_x'(v) \\ b_y'(u) = \lambda\,c_y'(v) \\ \end{cases} $$

with $0 < u < 1$, $0 < v < 2\pi$, $y_2 \in \mathbb{R}$, $\lambda \in \mathbb{R}$.

---

In light of this, defined the vector function $\mathbf{f} : \mathbb{R}^4 \to \mathbb{R}^4$ of law:

$$ \mathbf{f}(\mathbf{z}) := \left(b_x(u) - c_x(v), \; b_y(u) - c_y(v), \; b_x'(u) - \lambda\,c_x'(v), \; b_y'(u) - \lambda\,c_y'(v)\right) $$

with $\mathbf{z} \equiv \left(u,\,v,\,y_2,\,\lambda\right)$, the zeros of this function:

$$ \mathbf{f}(\mathbf{z}) = \mathbf{0} $$

can be calculated using the Newton-Raphson method, according to which:

$$ \mathbf{z}_{i + 1} = \mathbf{z}_i - \mathbf{J}_{\mathbf{f}}^{-1}(\mathbf{z}_i) \cdot \mathbf{f}(\mathbf{z}_i) $$

where $\mathbf{z}_0$ is an initial vector sufficiently close to the desired zero; the method stops when:

$$ \left\lVert \mathbf{J}_{\mathbf{f}}^{-1}(\mathbf{z}_i) \cdot \mathbf{f}(\mathbf{z}_i) \right\rVert \le \epsilon\,, $$

with $\epsilon > 0$ pre-established according to your needs.

---

Implementing this method (with some tricks) in Wolfram Mathematica 12.2:

{x0, y0} = {0, 10};
{x1, y1} = {0, 1.96};
x2 = 8;
{x3, y3} = {8, 0};
{xC, yC} = {6.713, 5.04};
r = 1.1;
ε = 10^-10;
countmax = 100;

bx[u_] = x0 (1 - u)^3 + 3 x1 u (1 - u)^2 + 3 x2 u^2 (1 - u) + x3 u^3;
by[u_] = y0 (1 - u)^3 + 3 y1 u (1 - u)^2 + 3 y2 u^2 (1 - u) + y3 u^3;

cx[v_] = xC + r Cos[v];
cy[v_] = yC + r Sin[v];

f[{u_, v_, y2_, λ_}] = {bx[u] - cx[v], 
                        by[u] - cy[v], 
                        bx'[u] - λ cx'[v], 
                        by'[u] - λ cy'[v]};

invJf[{u_, v_, y2_, λ_}] = Inverse[Grad[f[{u, v, y2, λ}], {u, v, y2, λ}]];

z = {1/2, 3π/2, (y0 + y1 + y3)/3, 10};
count = 0;
While[Norm[invJf[z].Transpose[{f[z]}]] > ε && count < countmax,
      z = z - Transpose[invJf[z].Transpose[{f[z]}]][[1]];
      If[z[[1]] < 0 || z[[1]] > 1, z[[1]] = 1/2];
      If[z[[2]] < π || z[[2]] > 2π, z[[2]] = 3π/2];
      count = count + 1
     ];

Print[{count, z}]

we get:

{6, {0.689356, 4.18946, 7.66843, 10.7859}}

that is, with only $6$ iterations, we obtain $y_2 = 7.66843$, which is what is desired.

Finally, to verify what has been obtained, by writing:

P0 = {x0, y0};
P1 = {x1, y1};
P2 = {x2, z[[3]]};
P3 = {x3, y3};

P0lab = P0 + {+0.3, +0.2};
P1lab = P1 + {+0.3, -0.2};
P2lab = P2 + {+0.3, +0.2};
P3lab = P3 + {+0.3, -0.2};

plot1 = ParametricPlot[Evaluate[{bx[u], by[u]} /. y2 -> z[[3]]], 
                       {u, 0, 1}, PlotStyle -> Red];

plot2 = ParametricPlot[{cx[v], cy[v]}, {v, 0, 2π}, PlotStyle -> Green];

plot3 = Graphics[{Directive[Blue, Dashed], Line[{P0, P1, P2, P3}]}];

Show[{plot1, plot2, plot3}, 
     AxesLabel -> {"x", "y"}, GridLines -> Automatic, PlotRange -> All,
     Epilog -> {Text[Style[Subscript["P", "0"], Blue], P0lab], 
                Text[Style[Subscript["P", "1"], Blue], P1lab],
                Text[Style[Subscript["P", "2"], Blue], P2lab], 
                Text[Style[Subscript["P", "3"], Blue], P3lab]}
    ]

we get:

which actually shows what you want. I hope it's clear enough, good luck! ^_^

---

---

Addendum 1

A fairly faithful approach to the algorithm implemented in Wolfram Mathematica 12.2 is:

1. define the parameters $x_0,\,y_0,\,x_1,\,y_1,\,x_2,\,x_3,\,y_3,\,x_C,\,y_C,\,r,\,\epsilon,\,\text{countmax}$;

2. initialize $\text{count} = 0$, $u = \dots$, $v = \dots$, $y_2 = \dots$, $\lambda = \dots$, then define $\mathbf{z} = \left(u,\,v,\,y_2,\,\lambda\right)^t$;

3. calculate the components: $$ \small \begin{aligned} & f_1 = x_0 + 3\left(x_1 - x_0\right)u + 3\left(x_2 - 2\,x_1 + x_0\right)u^2 + \left(x_3 - 3\,x_2 + 3\,x_1 - x_0\right)u^3 - x_C - r\,\cos v \\ & f_2 = y_0 + 3\left(y_1 - y_0\right)u + 3\left(y_2 - 2\,y_1 + y_0\right)u^2 + \left(y_3 - 3\,y_2 + 3\,y_1 - y_0\right)u^3 - y_C - r\,\sin v \\ & f_3 = 3\left(x_1 - x_0\right) + 6\left(x_2 - 2\,x_1 + x_0\right)u + 3\left(x_3 - 3\,x_2 + 3\,x_1 - x_0\right)u^2 + \lambda\,r\,\sin v \\ & f_4 = 3\left(y_1 - y_0\right) + 6\left(y_2 - 2\,y_1 + y_0\right)u + 3\left(y_3 - 3\,y_2 + 3\,y_1 - y_0\right)u^2 - \lambda\,r\,\cos v \\ \end{aligned} $$ then define $\mathbf{f} = \left(f_1,\,f_2,\,f_3,\,f_4\right)^t$;

4. calculate the components: $$ \begin{aligned} & J_{11} = 3\left(x_1 - x_0\right) + 6\left(x_2 - 2\,x_1 + x_0\right)u + 3\left(x_3 - 3\,x_2 + 3\,x_1 - x_0\right)u^2 \\ & J_{12} = r\,\sin v \\ & J_{13} = 0 \\ & J_{14} = 0 \\ \\ & J_{21} = 3\left(y_1 - y_0\right) + 6\left(y_2 - 2\,y_1 + y_0\right)u + 3\left(y_3 - 3\,y_2 + 3\,y_1 - y_0\right)u^2 \\ & J_{22} = -r\,\cos v \\ & J_{23} = 3\,u^2\,(1 - u) \\ & J_{24} = 0 \\ \\ & J_{31} = 6\left(x_2 - 2\,x_1 + x_0\right) + 6\left(x_3 - 3\,x_2 + 3\,x_1 - x_0\right)u \\ & J_{32} = \lambda\,r\,\cos v \\ & J_{33} = 0 \\ & J_{34} = r\,\sin v \\ \\ & J_{41} = 6\left(y_2 - 2\,y_1 + y_0\right) + 6\left(y_3 - 3\,y_2 + 3\,y_1 - y_0\right)u \\ & J_{42} = \lambda\,r\,\sin v \\ & J_{43} = 3\,u\,(2 - 3\,u) \\ & J_{44} = -r\,\cos v \\ \end{aligned} $$ then define $\mathbf{Jf} = \begin{pmatrix} J_{11} & J_{12} & J_{13} & J_{14} \\ J_{21} & J_{22} & J_{23} & J_{24} \\ J_{31} & J_{32} & J_{33} & J_{34} \\ J_{41} & J_{42} & J_{43} & J_{44} \\ \end{pmatrix}$;

5. calculate: $$ \begin{aligned} & \mathbf{invJf} = \text{inverse}(\mathbf{Jf})\,; \\ & \mathbf{resf} = \mathbf{invJf}.\mathbf{f} \; \text{[product row by column]}\,; \\ & \text{errf} = \sqrt{\mathbf{resf}_1^2 + \mathbf{resf}_2^2 + \mathbf{resf}_3^2 + \mathbf{resf}_4^2}\,; \end{aligned} $$

6. while $\text{errf} > \epsilon$ and $\text{count} < \text{countmax}$: $$ \begin{aligned} & \mathbf{z} = \mathbf{z} - \mathbf{resf}\,; \\ & \left(u,\,v,\,y_2,\,\lambda\right)^t = \mathbf{z}\,; \\ & \text{if} \; u < 0 \; \; \text{or} \; \; u > 1 \; \; \; \text{then} \; u = \frac{1}{2}\,; \\ & \text{if} \; v < \pi \; \; \text{or} \; \; v > 2\pi \; \; \text{then} \; v = \frac{3\pi}{2}\,; \\ & \text{repeat steps 3, 4, 5}\,; \\ & \text{count} = \text{count} + 1\,. \\ \end{aligned} $$

In particular, to make the calculations faster, it's possible to calculate the residual vector once and for all, reducing everything to a series of elementary operations that are easy to iterate:

(*Input data*)
{x0, y0} = {0, 10};
{x1, y1} = {0, 1.96};
x2 = 8;
{x3, y3} = {8, 0};
{xC, yC} = {6.713, 5.04};
r = 1.1;
ε = 10^-10;
countmax = 100;

(*Residual function definition*)
res[u_, v_, y2_, λ_] := Module[{f1, f2, f3, f4, J11, J12, J21, J22, J23,
         J31, J32, J34, J41, J42, J43, J44, num1, num2, num3, num4, den},
  
  f1 = x0 + 3 (x1 - x0) u + 3 (x2 - 2 x1 + x0) u^2 + (x3 - 3 x2 + 3 x1 - x0) u^3 - xC - r Cos[v];
  f2 = y0 + 3 (y1 - y0) u + 3 (y2 - 2 y1 + y0) u^2 + (y3 - 3 y2 + 3 y1 - y0) u^3 - yC - r Sin[v];
  f3 = 3 (x1 - x0) + 6 (x2 - 2 x1 + x0) u + 3 (x3 - 3 x2 + 3 x1 - x0) u^2 + λ r Sin[v];
  f4 = 3 (y1 - y0) + 6 (y2 - 2 y1 + y0) u + 3 (y3 - 3 y2 + 3 y1 - y0) u^2 - λ r Cos[v];
  
  J11 = 3 (x1 - x0) + 6 (x2 - 2 x1 + x0) u + 3 (x3 - 3 x2 + 3 x1 - x0) u^2;
  J12 = r Sin[v];

  J21 = 3 (y1 - y0) + 6 (y2 - 2 y1 + y0) u + 3 (y3 - 3 y2 + 3 y1 - y0) u^2;
  J22 = -r Cos[v];
  J23 = 3 u^2 (1 - u);

  J31 = 6 (x2 - 2 x1 + x0) + 6 (x3 - 3 x2 + 3 x1 - x0) u;
  J32 = λ r Cos[v];
  J34 = r Sin[v];

  J41 = 6 (y2 - 2 y1 + y0) + 6 (y3 - 3 y2 + 3 y1 - y0) u;
  J42 = λ r Sin[v];
  J43 = 3 u (2 - 3 u);
  J44 = -r Cos[v];
  
  num1 = f1 (J23 (J34 J42 - J32 J44) - J22 J34 J43) +
         J12 (f2 J34 J43 + f3 J23 J44 - f4 J23 J34);
  
  num2 = f1 (J23 (J31 J44 - J34 J41) + J21 J34 J43) -
         J11 (f2 J34 J43 + f3 J23 J44 - f4 J23 J34);
  
  num3 = (f1 J21 - f2 J11) (J32 J44 - J34 J42) -
         (f1 J22 - f2 J12) (J31 J44 - J34 J41) + 
         (f3 J44 - f4 J34) (J11 J22 - J12 J21);
  
  num4 = (f1 J31 - f3 J11) (J22 J43 - J23 J42) -
         (f1 J32 - f3 J12) (J21 J43 - J23 J41) +
         (f2 J43 - f4 J23) (J11 J32 - J12 J31);
  
  den = J11 (J23 (J34 J42 - J32 J44) - J22 J34 J43) +
        J12 (J23 (J31 J44 - J34 J41) + J21 J34 J43);
  
  Return[{num1, num2, num3, num4} / den]];

(*Attempt solution*)
{u, v, y2, λ} = {1/2, 3π/2, (y0 + y1 + y3)/3, 10};
vres = res[u, v, y2, λ];
errq = vres[[1]]^2 + vres[[2]]^2 +
       vres[[3]]^2 + vres[[4]]^2;
count = 0;

(*Update solution*)
While[errq > ε^2 && count < countmax,
     {u, v, y2, λ} = {u, v, y2, λ} - vres;
     If[u < 0 || u > 1, u = 1/2];
     If[v < π || v > 2π, v = 3π/2];
     vres = res[u, v, y2, λ];
     errq = vres[[1]]^2 + vres[[2]]^2 +
            vres[[3]]^2 + vres[[4]]^2;
     count = count + 1
    ];

(*Output data*)
Print[{count, {u, v, y2, λ}}]

Although the definition of these parameters is very tedious, once this work is done, there is the advantage of being able to dominate the algorithm in an optimal way, as there are no more hidden steps.

In particular, I would like to underline the fact that the high efficiency of the Newton-Raphson method is paid for in terms of low robustness, since the convergence of this method is highly influenced by the choice of the first attempt solution vector.

Precisely for this reason, based on the problem under examination, it's necessary to intervene by optimizing it through some more or less intelligent measures. In this case, I intervened inside the while loop with a couple of if conditional checks to make $u$ and $v$ converge to a value within the respective intervals, but this can certainly be done more efficiently, as well as you will almost certainly have to intervene in the testing to adapt it to your needs.

I hope I have explained myself sufficiently, good job! ^_^

---

---

Addendum 2

If the specific problem doesn't allow to identify a robust strategy in the choice of the initial vector to trigger the Newton-Raphson method, it's necessary to change direction.

In particular, it's possible to abandon the world of local convergence methods by migrating to that of global convergence methods in which, when it comes to polynomials, the Aberth-Ehrlich method is undoubtedly an excellent compromise between algorithmic simplicity and at the same time efficiency/robustness.

On an algorithmic level, I would proceed like this:

1. through simple substitutions, the above system of equations can be simplified as follows: $$ \begin{cases} a + b\,u + c\,u^2 + d\,u^3 = \cos v \\ e + f\,u + \left(g + 3\,y_2/r\right)u^2 + \left(h - 3\,y_2/r\right)u^3 = \sin v \\ b + 2\,c\,u + 3\,d\,u^2 = -\lambda\,\sin v \\ f + 2\left(g + 3\,y_2/r\right)u + 3\left(h - 3\,y_2/r\right)u^2 = \lambda\,\cos v \\ \end{cases} $$ with $0 < u < 1$, $0 < v < 2\pi$, $y_2 \in \mathbb{R}$, $\lambda \in \mathbb{R}$;

2. dividing the first two equations member to member, I obtain the expression of $\sin v/\cos v$, then dividing the last two equations I obtain the expression of $-\sin v/\cos v$ which, when suitably equalized, allow to obtain a quadratic equation in $y_2$, an unknown which can therefore be expressed as a function of $u$;

3. adding the squares member to member of the first two equations, substituting the expression of $y_2$ and simplifying everything with further squares, it's possible to obtain a polynomial equation $p(u) = 0$, which in the most general case is 12th degree, where the coefficient of the monomial of maximum degree can be made unitary by dividing all the others by itself;

4. the time has therefore come to bring heavy artillery into the field by applying the above numerical method, which is conceived precisely for monical polynomial equations, in general with complex coefficients, through which it's possible to calculate all twelve complex roots simultaneously by updating them with the following formulation: $$ u_j' = u_j - \frac{1}{\frac{p'\left(u_j\right)}{p\left(u_j\right)} - \begin{aligned}\sum_{k = 1, \, k \ne j}^{12}\end{aligned} \frac{1}{u_j - u_k}} \;; $$

5. all that remains is to examine these roots one by one by selecting only those with a real part between $0$, $1$ and with an imaginary part close to zero, then through simple checks it's possible to determine the set of values $\left(u,\,v,\,y_2,\,\lambda\right)$ that satisfy all four equations placed in the system (due to the way the problem has been formulated, there may be more solutions than there may be, unless we introduce further constraints on the input parameters).

The advantage of this other method is that it doesn't require a particularly reasoned initial vector, it's sufficient to define it once and for all so that the (complex) components are asymmetrical and uniformly distributed (generally assumed on a circle).

On the other hand, it's necessary to refer to the arithmetic of complex numbers (which is the lesser evil, since it's easily algorithmizable), so it's necessary to establish a tolerance with which to discriminate the real roots from the complex roots by comparing it with the imaginary part (which could cause problems, since it isn't so obvious to define it arbitrarily).

Anyway, I tried to write this other algorithm in Visual Basic for Applications (VBA) in Excel and I was very satisfied, it seems very robust, although obviously it can always be improved:

Function sum(x, y)
sum = Array(x(0) + y(0), x(1) + y(1))
End Function

Function diff(x, y)
diff = Array(x(0) - y(0), x(1) - y(1))
End Function

Function prod(x, y)
prod = Array(x(0) * y(0) - x(1) * y(1), x(1) * y(0) + x(0) * y(1))
End Function

Function div(x, y)
div = Array((x(0) * y(0) + x(1) * y(1)) / (y(0) ^ 2 + y(1) ^ 2), (x(1) * y(0) - x(0) * y(1)) / (y(0) ^ 2 + y(1) ^ 2))
End Function

Function pow(x, n)
pow = Array(Cos(n * WorksheetFunction.Atan2(x(0), x(1))) * (x(0) ^ 2 + x(1) ^ 2) ^ (n / 2), Sin(n * WorksheetFunction.Atan2(x(0), x(1))) * (x(0) ^ 2 + x(1) ^ 2) ^ (n / 2))
End Function

Sub solve()

'IMPORT DATA'
x0 = Cells(4, 4)
y0 = Cells(5, 4)
x1 = Cells(7, 4)
y1 = Cells(8, 4)
x2 = Cells(10, 4)
x3 = Cells(13, 4)
y3 = Cells(14, 4)
xC = Cells(16, 4)
yC = Cells(17, 4)
r = Cells(18, 4)
imax = Cells(26, 4)
tol = Cells(27, 4)
segno = Cells(28, 4)

'SCALING'
a = (x0 - xC) / r
b = 3 * (x1 - x0) / r
c = 3 * (x2 - 2 * x1 + x0) / r
d = (x3 - 3 * x2 + 3 * x1 - x0) / r
e = (y0 - yC) / r
f = 3 * (y1 - y0) / r
g = 3 * (-2 * y1 + y0) / r
h = (y3 + 3 * y1 - y0) / r

'POLYNOMIAL COEFFICIENTS'
u0 = 4 * (-1 + a ^ 2) * (-1 + a ^ 2 + e ^ 2)
u1 = -4 * (3 + 3 * a ^ 4 - 3 * a ^ 3 * b + a * b * (3 - 2 * e ^ 2) + e * (-3 * e + f) + a ^ 2 * (-6 + 3 * e ^ 2 - e * f))
u2 = 9 * a ^ 4 + a ^ 3 * (-38 * b + 8 * c) + (-9 + 4 * b ^ 2) * (-1 + e ^ 2) + 8 * a * c * (-1 + e ^ 2) + 14 * e * f - f ^ 2 + a ^ 2 * (-18 + 13 * b ^ 2 + 9 * e ^ 2 - 14 * e * f + f ^ 2) + a * b * (38 + 8 * e * (-3 * e + f))
u3 = 2 * (a ^ 3 * (15 * b + 2 * (-7 * c + d)) + 2 * b * c * (-1 + 2 * e ^ 2) + b ^ 2 * (7 + 2 * e * (-3 * e + f)) + a * (3 * b ^ 3 + 14 * c - 2 * d + b * (-15 + 9 * e ^ 2 - 14 * e * f + f ^ 2) + 4 * e * (d * e + c * (-3 * e + f))) + 2 * (f ^ 2 + e * (-3 * f + g + h)) - 2 * a ^ 2 * (b * (11 * b - 4 * c) + f ^ 2 + e * (-3 * f + g + h)))
u4 = -22 * a * b ^ 3 + b ^ 4 + b ^ 2 * (-12 + 37 * a ^ 2 + 10 * a * c + 9 * e ^ 2 - 14 * e * f + f ^ 2) + 2 * (3 * a ^ 3 * (4 * c - 3 * d) + 2 * c ^ 2 * e ^ 2 + a * (c * (-12 + 9 * e ^ 2 - 14 * e * f + f ^ 2) + d * (9 + 4 * e * (-3 * e + f))) - 3 * e * (g + h) + f * (-2 * f + g + h) + a ^ 2 * (2 * c ^ 2 + 2 * f ^ 2 + 3 * e * (g + h) - f * (g + h))) + b * (6 * a ^ 2 * d + 8 * d * e ^ 2 + c * (18 - 62 * a ^ 2 + 8 * e * (-3 * e + f)) - 8 * a * (f ^ 2 + e * (-3 * f + g + h)))
u5 = 2 * (-2 * b ^ 4 + b ^ 3 * c + 9 * a ^ 3 * d + b * (c * (-9 + 9 * e ^ 2 - 14 * e * f + f ^ 2) + 4 * d * (1 + e * (-3 * e + f))) + 2 * (c * (d + 2 * d * e ^ 2) + c ^ 2 * (1 + e * (-3 * e + f)) - f * (g + h)) + a ^ 2 * (29 * b * c - 10 * c ^ 2 - 18 * b * d + 2 * f * (g + h)) - 2 * b ^ 2 * (f ^ 2 + e * (-3 * f + g + h)) + a * (10 * b ^ 3 + b ^ 2 * (-22 * c + d) + d * (-9 + 9 * e ^ 2 - 14 * e * f + f ^ 2) + 2 * b * (c ^ 2 + 2 * f ^ 2 + 3 * e * (g + h) - f * (g + h)) - 4 * c * (f ^ 2 + e * (-3 * f + g + h))))
u6 = 4 * b ^ 4 - 10 * b ^ 3 * c + 4 * d ^ 2 - 3 * a ^ 2 * d ^ 2 + 4 * d ^ 2 * e ^ 2 + 24 * a * d * e * f - 8 * a * d * f ^ 2 + c ^ 2 * (-6 + 22 * a ^ 2 + 9 * e ^ 2 - 14 * e * f + f ^ 2) - 8 * a * d * e * g - g ^ 2 + a ^ 2 * g ^ 2 - 8 * a * d * e * h - 2 * g * h + 2 * a ^ 2 * g * h - h ^ 2 + a ^ 2 * h ^ 2 + b ^ 2 * (46 * a * c + c ^ 2 - 22 * a * d + 4 * f ^ 2 + 6 * e * (g + h) - 2 * f * (g + h)) - 2 * c * (d * (1 + 9 * a ^ 2 + 12 * e ^ 2 - 4 * e * f) - 6 * a * e * (g + h) + 2 * a * f * (-2 * f + g + h)) + 2 * b * (21 * a ^ 2 * d + d * (-6 + 9 * e ^ 2 - 14 * e * f + f ^ 2) + a * (-13 * c ^ 2 - 2 * c * d + 4 * f * (g + h)) - 4 * c * (f ^ 2 + e * (-3 * f + g + h)))
u7 = 2 * (b ^ 3 * (6 * c - 2 * d) - 3 * c * d + 15 * a ^ 2 * c * d - 3 * d ^ 2 + 9 * c * d * e ^ 2 - 6 * d ^ 2 * e ^ 2 + 6 * c ^ 2 * e * f - 14 * c * d * e * f + 2 * d ^ 2 * e * f - 2 * c ^ 2 * f ^ 2 + c * d * f ^ 2 - 2 * c ^ 2 * e * g - 2 * c ^ 2 * e * h + b ^ 2 * (16 * a * d - c * (4 * c + d) + 2 * f * (g + h)) - 2 * a * (c ^ 3 + c ^ 2 * d - 3 * d * e * (g + h) - 2 * c * f * (g + h) + d * f * (-2 * f + g + h)) + b * (2 * c * (2 * f ^ 2 + 3 * e * (g + h) - f * (g + h)) + a * (17 * c ^ 2 - 8 * c * d - 3 * d ^ 2 + (g + h) ^ 2) - 4 * d * (f ^ 2 + e * (-3 * f + g + h))))
u8 = 8 * b ^ 3 * d + 9 * a ^ 2 * d ^ 2 + 9 * d ^ 2 * e ^ 2 + 24 * c * d * e * f - 14 * d ^ 2 * e * f + 4 * c ^ 2 * f ^ 2 - 8 * c * d * f ^ 2 + d ^ 2 * f ^ 2 + 6 * c ^ 2 * e * g - 8 * c * d * e * g - 2 * c ^ 2 * f * g - 2 * c * (4 * d * e + c * (-3 * e + f)) * h + b ^ 2 * (13 * c ^ 2 - 2 * c * d - 2 * d ^ 2 + (g + h) ^ 2) + 2 * b * (-c ^ 3 - c ^ 2 * d + d * (3 * a * d + 4 * f ^ 2 + 6 * e * (g + h) - 2 * f * (g + h)) + c * (22 * a * d + 4 * f * (g + h))) + 2 * a * (4 * c ^ 3 + c ^ 2 * d + 4 * d * f * (g + h) + c * (-3 * d ^ 2 + (g + h) ^ 2))
u9 = 2 * (2 * b ^ 2 * d * (4 * c + d) + 2 * c ^ 2 * f * (g + h) + 2 * c * d * (2 * f ^ 2 + 3 * e * (g + h) - f * (g + h)) + a * d * ((7 * c - d) * (c + d) + (g + h) ^ 2) + b * (3 * c ^ 3 + 2 * c ^ 2 * d + 6 * a * d ^ 2 - c * d ^ 2 + 4 * d * f * (g + h) + c * (g + h) ^ 2) - 2 * d ^ 2 * (f ^ 2 + e * (-3 * f + g + h)))
u10 = c ^ 4 + 2 * c ^ 3 * d + 2 * c * d * (3 * a * d + 5 * b * d + 4 * f * (g + h)) + c ^ 2 * (10 * b * d + d ^ 2 + (g + h) ^ 2) + 2 * d * (2 * b ^ 2 * d + b * (g + h) ^ 2 + d * (3 * a * d + 2 * f ^ 2 + 3 * e * (g + h) - f * (g + h)))
u11 = 2 * d * (c ^ 3 + 2 * c ^ 2 * d + 2 * d * (b * d + f * (g + h)) + c * (2 * b * d + d ^ 2 + (g + h) ^ 2))
u12 = d ^ 2 * ((c + d) ^ 2 + (g + h) ^ 2)

'VECTORING'
ReDim coeff(1 To 1, 1 To 13)
coeff(1, 1) = Array(u0, 0)
coeff(1, 2) = Array(u1, 0)
coeff(1, 3) = Array(u2, 0)
coeff(1, 4) = Array(u3, 0)
coeff(1, 5) = Array(u4, 0)
coeff(1, 6) = Array(u5, 0)
coeff(1, 7) = Array(u6, 0)
coeff(1, 8) = Array(u7, 0)
coeff(1, 9) = Array(u8, 0)
coeff(1, 10) = Array(u9, 0)
coeff(1, 11) = Array(u10, 0)
coeff(1, 12) = Array(u11, 0)
coeff(1, 13) = Array(u12, 0)

'COEFFICIENTS ANALYSIS'
n = 13
flag = 0
While flag = 0
    If Abs(coeff(1, n)(0)) < 10 ^ -15 Then
        n = n - 1
    Else
        flag = 1
    End If
Wend

'NORMALIZE COEFFICIENTS'
For i = 1 To n - 1
    coeff(1, i) = div(coeff(1, i), coeff(1, n))
Next i
n = n - 1

'ATTEMPT u-SOLUTION'
ReDim sol(1 To 1, 1 To n)
For i = 1 To n
    sol(1, i) = Array(Cos(2 * WorksheetFunction.Pi() * i / n), Sin(2 * WorksheetFunction.Pi() * i / n))
Next i

'UPDATE u-SOLUTION'
For i = 1 To imax
    
    'ABERTH-EHRLICH METHOD'
    For j = 1 To n
        
        'GENERATE p'
        p = Array(0, 0)
        For k = 1 To n
            p = sum(p, prod(coeff(1, k), pow(sol(1, j), k - 1)))
        Next k
        p = sum(p, pow(sol(1, j), n))
        
        'GENERATE dp'
        dp = Array(0, 0)
        For k = 2 To n
            dp = sum(dp, prod(Array(k - 1, 0), prod(coeff(1, k), pow(sol(1, j), k - 2))))
        Next k
        dp = sum(dp, prod(Array(n, 0), pow(sol(1, j), n - 1)))
        
        'GENERATE q'
        q = Array(0, 0)
        For k = 1 To n
            If k <> j Then
                q = sum(q, div(Array(1, 0), diff(sol(1, j), sol(1, k))))
            End If
        Next k
        
        sol(1, j) = diff(sol(1, j), div(p, diff(dp, prod(p, q))))

    Next j

Next i

'SEARCH y2-SOLUTION'
i = 1
flag = 0
While flag = 0 And i <= n

    If sol(1, i)(0) > 0 And sol(1, i)(0) < 1 And Abs(sol(1, i)(1)) < tol Then
        
        u = sol(1, i)(0)
        a0 = 9 * u ^ 3 * (2 - 3 * u) * (1 - u)
        b0 = 3 * u * ((2 - 3 * u) * (e + f * u + g * u ^ 2 + h * u ^ 3) + u * (1 - u) * (f + 2 * g * u + 3 * h * u ^ 2))
        c0 = (a + b * u + c * u ^ 2 + d * u ^ 3) * (b + 2 * c * u + 3 * d * u ^ 2) + (e + f * u + g * u ^ 2 + h * u ^ 3) * (f + 2 * g * u + 3 * h * u ^ 2)
        cosine = a + b * u + c * u ^ 2 + d * u ^ 3

        If b0 ^ 2 - 4 * a0 * c0 > 0 And Abs(2 * a0) > 10 ^ -15 And cosine > -1 And cosine < 1 Then
            
            y2 = r * (-b0 - (b0 ^ 2 - 4 * a0 * c0) ^ (1 / 2)) / (2 * a0)
            sine = e + f * u + (g + 3 * y2 / r) * u ^ 2 + (h - 3 * y2 / r) * u ^ 3
            v = (1 - segno) * WorksheetFunction.Pi() + segno * WorksheetFunction.Acos(cosine)
            
            If Abs(sine - Sin(v)) < tol Then
                flag = 1
                lambda = -(b + 2 * c * u + 3 * d * u ^ 2) / sine
            End If
                
            If flag = 0 Then

                y2 = r * (-b0 + (b0 ^ 2 - 4 * a0 * c0) ^ (1 / 2)) / (2 * a0)
                sine = e + f * u + (g + 3 * y2 / r) * u ^ 2 + (h - 3 * y2 / r) * u ^ 3
            
                If Abs(sine - Sin(v)) < tol Then
                    flag = 1
                    lambda = -(b + 2 * c * u + 3 * d * u ^ 2) / sine
                End If
                
            End If
            
        End If
        
    End If
    
    i = i + 1
    
Wend

'PRINT SOLUTION'
If flag = 1 Then
    Cells(4, 21) = u
    Cells(5, 21) = v
    Cells(11, 4) = y2
    Cells(6, 21) = y2
    Cells(7, 21) = lambda
Else
    Cells(4, 21) = ""
    Cells(5, 21) = ""
    Cells(11, 4) = ""
    Cells(6, 21) = ""
    Cells(7, 21) = ""
    Output = MsgBox("Solution not found!", vbCritical, "WARNING!")
End If

'OPTIMIZE ASPECT RATIO'
xmin = Application.WorksheetFunction.Min(Union(Range("D4:D4"), Range("D7:D7"), Range("D10:D10"), Range("D13:D13"), Range("F4:F104"), Range("I4:I104")))
xmax = Application.WorksheetFunction.Max(Union(Range("D4:D4"), Range("D7:D7"), Range("D10:D10"), Range("D13:D13"), Range("F4:F104"), Range("I4:I104")))
ymin = Application.WorksheetFunction.Min(Union(Range("D5:D5"), Range("D8:D8"), Range("D11:D11"), Range("D14:D14"), Range("G4:G104"), Range("J4:J104")))
ymax = Application.WorksheetFunction.Max(Union(Range("D5:D5"), Range("D8:D8"), Range("D11:D11"), Range("D14:D14"), Range("G4:G104"), Range("J4:J104")))

minabs = Application.WorksheetFunction.Min(xmin, ymin) * 110 / 100
maxabs = Application.WorksheetFunction.Max(xmax, ymax) * 110 / 100

ActiveSheet.ChartObjects("Grafico 1").Activate
ratio = 5 / 6

If xmax - xmin > ratio * (ymax - ymin) Then
    ActiveChart.Axes(xlCategory).MinimumScale = minabs
    ActiveChart.Axes(xlCategory).MaximumScale = maxabs
    ActiveChart.Axes(xlValue).MinimumScale = minabs / ratio
    ActiveChart.Axes(xlValue).MaximumScale = maxabs / ratio
Else
    ActiveChart.Axes(xlCategory).MinimumScale = minabs * ratio
    ActiveChart.Axes(xlCategory).MaximumScale = maxabs * ratio
    ActiveChart.Axes(xlValue).MinimumScale = minabs
    ActiveChart.Axes(xlValue).MaximumScale = maxabs
End If

ActiveSheet.Range("A1").Select

End Sub

I hope I haven't written nonsense, in case let me know in the comments. Good luck. ^_^

Answer 2

Let us plug the parametric equation of the Bezier into the implicit equation of the circle:

$$\left(\sum_{k=0}^3 x_kB_k(t)-a\right)^2+\left(\sum_{k=0}^3 y_kB_k(t)-b\right)^2=r^2$$ where the $B_k$'s denote the Bernstein polynomials.

We express that the curve is tangent, or that the above equation has a double root in $t$:

$$\sum_{k=0}^3 x_kB'_k(t)\left(\sum_{k=0}^3 x_kB_k(t)-a\right)+\sum_{k=0}^3 y_kB'_k(t)\left(\sum_{k=0}^3 y_kB_k(t)-b\right)=0.$$

This forms a system of two equations in two unknowns, namely $t$ and $y_2$.

If you expand, you see that the two equations are quadratic in $y_2$. Then writing the condition for two quadratic equations to share a root, which is quartic in the coefficients, you can eliminate $y_2$.

In the end, you will obtain a monstrous polynomial in $t$ that could be of a degree up to $60$ (needs to be checked). This is a clear indication that your problem is arduous and must be solved numerically.

---

In a more compact form, the equations are

$$(Y(t)+B_2(t)y_2)^2+X^2(t)-r^2=0$$ $$B'_2(t)y_2(Y(t)+B_2(t)y_2)-X'(t)X(t)=0$$ where $X,Y$ denote cubic polynomials.