# Cubic Bezier curve tangent to Circle

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

• Are $(x_0,y_0)$ and $(x_3,y_3)$ given fixed points, $x_0=x_1$, $x_2=x_3$, and the question is to find corresponding $y_2$ for any $y_1$? Also, it looks like wlog you can simplify the question by setting $a=b=0,\ R=1$. Apr 27, 2021 at 14:02
• Yes, P0 and P3 are fixed points, and for simplicity I can say that x0 = x1 = y3 = 0. But what I don't understand is that if a = b = 0, then my circle is no longer in the right place, and its radius is not necessarily equal to 1 ... Apr 28, 2021 at 17:26
• Yes, you need to check which option results in simpler expressions.. You get $a = b = 0,\ R=1$ by subtracting $(a,b)$ from all the other coordinates and then scaling them all by $1/R$. Apr 28, 2021 at 17:38
• Ok ... Thanks.. Apr 29, 2021 at 8:17

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! ^_^

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! ^_^

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) = ""
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. ^_^

• Thank you very much for your very comprehensive response. Your implementation in Wolfram interests me a lot, because that's precisely what I want to do… but in Javascript ! I want to use this calculation to create shapes in SVG format (this format uses Bezier cubic curves), and all of this written in Javascript. In your code, I see things that look like JavaScript like variables (x2 = 8; or count = 0;), or the function "while()" that will surely determine the Newton-Raphson iterations. Unfortunately I don't know how to translate this code into Javascript. Do you have an idea ? Apr 28, 2021 at 17:34
• As you can see my math level is not very high…:-) I notably spent a lot of time looking for what this "$J_{f}^{-1}$" corresponded to in your answer, before discovering the "Jacobian Matrix" that I did not know… :-) I'm learning things… It's great ! Your addendum is really great and will help me a lot to finalize my problem. I feel that the solution is not far away ! Many thanks again for your time on this issue. Apr 30, 2021 at 16:19

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.

• Indeed, I had started by "ploting" the Bezier curve into the equation of the circle, and I quickly realized that I arrived at an algebraically insoluble sixth degree equation. I agree, the solution will come with a numerical method. Your compact form of equations interests me. I'll see where it takes me. Thanks. Apr 30, 2021 at 16:23