Closest points on two parabolas in 3D

Question

Two parabolas are given in parametric form

$r_1 = r_{01} + v_1 \ t_1 + v_2 \ t_1^2 $

$r_2 = r_{02} + w_1 \ t_2 + w_2 \ t_2^2 $

where $r_{01}, v_1, v_2, r_{02}, w_1, w_2 \in \mathbb{R}^3$

It is required to find $t_1$ and $t_2$ that will result in $r_1$ and $r_2$ being as close as possible.

My attempt

The squared distance function is

$ f = (r_1 - r_2) \cdot (r_1 - r_2) $

Differentiating with respect to $t_1$ and $t_2$, and setting them to $0$, gives the following two equations

$ (r_1 - r_2) \cdot (v_1 + 2 t_1 v_2 ) = 0 $

$ (r_1 - r_2) \cdot (w_1 + 2 t_2 w_2) = 0 $

These two equations involve $t_1, \ t_1^2, \ t_2, \ t_2^2, \ t_1^2 t_2 , \ t_1 t_2^2, \ t_1^3, \ t_2^3 $

So probably the best way to solve them for $t_1$ and $t_2$ is numerically, through the Newton-Raphson method. That is exactly what I did, and I got the solution in $5$ iterations which is typical for nonlinear systems of equations with $2$ variables.

The code I used is listed in my solution below.

Any comments, hints, or alternative solutions are welcome and highly appreciated.

Answer 1

Below is the code that I used to implement the Newton-Raphson iteration. It is written in Microsoft Excel VBA (Visual Basic for Applications) script.

Public Sub closest_points_on_two_parabolas_3D()

'(r01 - r02 + t1 v1 + t1^2 v2 - t2 w1 - t2^2 w2 ) . (v1 + 2 t1 v2) = 0

'(r01 - r02 + t1 v1 + t1^2 v2 - t2 w1 - t2^2 w2 ) . (w1 + 2 t2 w2) = 0


Dim v1(3), v2(3), w1(3), w2(3), r01(3), r02(3) As Double

Dim x(2), y(2), jac(2, 3) As Double

r01(1) = 0
r01(2) = 5
r01(3) = 6

r02(1) = 7
r02(2) = 10
r02(3) = 1

For i = 1 To 3
v1(i) = ActiveSheet.Cells(i, 1)
v2(i) = ActiveSheet.Cells(i, 2)
w1(i) = ActiveSheet.Cells(i, 3)
w2(i) = ActiveSheet.Cells(i, 4)
Next i

x(1) = 0
x(2) = 0
success = 0

For ic = 1 To 200

Call find_y_two_parabolas(r01, v1, v2, r02, w1, w2, x, y)


If norm(y, 2) < 0.0000001 Then

success = 1
MsgBox ("Converged in " + Str(ic - 1))
Exit For

End If

Call find_jac_two_parabolas(r01, v1, v2, r02, w1, w2, x, jac)

For i = 1 To 2
jac(i, 3) = -y(i)
Next i

Call reduce_matrix(jac, 2, 3, 0.00000001, det)

For i = 1 To 2
x(i) = x(i) + jac(i, 3)
Next i

Next ic

If success = 0 Then
MsgBox ("did not converge")
Exit Sub
End If

For i = 1 To 3
u1(i) = r01(i) + x(1) * v1(i) + x(1) ^ 2 * v2(i)
u2(i) = r02(i) + x(2) * w1(i) + x(2) ^ 2 * w2(i)
u3(i) = u2(i) - u1(i)
ActiveSheet.Cells(i + 5, 1) = u1(i)
ActiveSheet.Cells(i + 5, 2) = u2(i)
Next i

ActiveSheet.Cells(10, 1) = norm(u3, 2)

End Sub

---

Public Sub find_y_two_parabolas(r01, v1, v2, r02, w1, w2, x, y)
Dim r1(3), r2(3) As Double
Dim tv1(3), tv2(3) As Double

'(r01 - r02 + t1 v1 + t1^2 v2 - t2 w1 - t2^2 w2 ) . (v1 + 2 t1 v2) = 0

'(r01 - r02 + t1 v1 + t1^2 v2 - t2 w1 - t2^2 w2 ) . (w1 + 2 t2 w2) = 0

For i = 1 To 3
r1(i) = r01(i) + x(1) * v1(i) + x(1) ^ 2 * v2(i)
r2(i) = r02(i) + x(2) * w1(i) + x(2) ^ 2 * w2(i)
u1(i) = r1(i) - r2(i)

tv1(i) = v1(i) + 2 * x(1) * v2(i)
tv2(i) = w1(i) + 2 * x(2) * w2(i)
Next i

y(1) = dot(u1, tv1)
y(2) = dot(u1, tv2)


End Sub

---

Public Sub find_jac_two_parabolas(r01, v1, v2, r02, w1, w2, x, jac)
Dim r1(3), r2(3) As Double
Dim tv1(3), tv2(3) As Double

'(r01 - r02 + t1 v1 + t1^2 v2 - t2 w1 - t2^2 w2 ) . (v1 + 2 t1 v2) = 0

'(r01 - r02 + t1 v1 + t1^2 v2 - t2 w1 - t2^2 w2 ) . (w1 + 2 t2 w2) = 0

For i = 1 To 3
r1(i) = r01(i) + x(1) * v1(i) + x(1) ^ 2 * v2(i)
r2(i) = r02(i) + x(2) * w1(i) + x(2) ^ 2 * w2(i)
u1(i) = r1(i) - r2(i)

tv1(i) = v1(i) + 2 * x(1) * v2(i)
tv2(i) = w1(i) + 2 * x(2) * w2(i)
Next i

jac(1, 1) = dot(tv1, tv1) + 2 * dot(u1, v2)
jac(1, 2) = -dot(tv2, tv1)

jac(2, 1) = dot(tv1, tv2)
jac(2, 2) = -dot(tv2, tv2) + 2 * dot(w2, u1)

End Sub

---