TransWikia.com

Finding the intersection(s) of a 3D parametric surface and a line

Mathematica Asked by GalZoidberg on November 27, 2020

I have some 3D parametric curve par and some line line.

I want to find parameters u,v and t for the intersection points of par and line.

I tried to do that with NSolve in many different ways and always got the following error message:

"This system cannot be solved with the methods available to NSolve."

Below is just an example, but I’m looking for a method, that works for every parametric curve and every line.

Also I don’t want to get x,y and z but u,v and t

How to fix it?

Clear[u, v, t];
par = 5 {Cos[u], Cos[v] + Sin[u], Sin[v]};
line = {1, 2, 3} + {4, -5, 6} t;
NSolve[
 par == line
  && 0 < u < 2 [Pi]
  && -[Pi] < v < [Pi],
 {t, u, v}, Reals]

2 Answers

In general, Reduce is the most powerful Mathematica function for solving equations.

Reduce[par == line && 0 < u < 2 Pi && -Pi < v < Pi, {t, u, v}, Reals]
(* (t == AlgebraicNumber[Root[-18291750000 + 235480000*#1 - 237800*#1^2 + 
        48*#1^3 + #1^4 & , 1, 0], {0, 1/725, 0, 0}] && 
    u == 2*ArcTan[AlgebraicNumber[Root[-18291750000 + 235480000*#1 - 
        237800*#1^2 + 48*#1^3 + #1^4 & , 1, 0], 
        {3961/4536, -2479/1315440, 1121/353220000, 83/19073880000}]] && 
    v == 2*ArcTan[AlgebraicNumber[Root[-18291750000 + 235480000*#1 - 
        237800*#1^2 + 48*#1^3 + #1^4 & , 1, 0], {3505/1484, -9221/6455400, 
        17053/9360330000, 197/18720660000}]]) || 
    (t == AlgebraicNumber[Root[-18291750000 + 235480000*#1 - 237800*#1^2 + 
        48*#1^3 + #1^4 & , 2, 0], {0, 1/725, 0, 0}] && 
     u == 2*ArcTan[AlgebraicNumber[Root[-18291750000 + 235480000*#1 - 
        237800*#1^2 + 48*#1^3 + #1^4 & , 2, 0], 
        {3961/4536, -2479/1315440, 1121/353220000, 83/19073880000}]] && 
     v == 2*ArcTan[AlgebraicNumber[Root[-18291750000 + 235480000*#1 - 
         237800*#1^2 + 48*#1^3 + #1^4 & , 2, 0], 
         {3505/1484, -9221/6455400, 17053/9360330000, 197/18720660000}]]) @)

% // N
(* (t == -1.07503 && u == 2.29164 && v == -0.761533) || 
   (t == 0.116633 && u == 1.27311 && v == 2.30858) *)

The intersections can be visualized by

Show[
    ParametricPlot3D[par, {u, 0, 2 Pi}, {v, -Pi, Pi}, 
        PlotStyle -> Opacity[.5], LabelStyle -> {15, Bold, Black}],
    ParametricPlot3D[line, {t, -2, 1}, PlotStyle -> {Black, Thick}],
    ListPointPlot3D[{line /. t -> -1.07503, line /. t -> 0.11663}, PlotStyle -> Red]]

enter image description here

Addendum: Use FindRoot

In the event that Reduce does not provide a solution, FindRoot almost always will, but requires multiple initial guesses to obtain multiple intersections, as is the case here.

FindRoot[par == line, {t, 0}, {u, Pi}, {v, 0}]
FindRoot[par == line, {t, 0}, {u, Pi}, {v, 2}]
(* {t -> -1.07503, u -> 2.29164, v -> -0.761533} *)
(* {t -> 0.116633, u -> 1.27311, v -> 2.30858} *)

Answered by bbgodfrey on November 27, 2020

In general, first try Reduce as already mentioned. But some more difficult curves will not reduce or Mathematica will hang. For example:

par = {1, Erf[t], Cos[t]};
line = {1, 1/2, 0} + {0, 1, 1} s;
Reduce[par == line && 0 < s < 2 && 0 < t < 2, {s, t}]

Instead you can try minimizing the distance between the line and the curve:

result = NMinimize[{SquaredEuclideanDistance[par, line], 0 < s < 2, 0 < t < 2}, {s, t}]
(* result: {4.42622*10^-19, {s -> 0.399155, t -> 1.1602}} *)
{par,line} /. Last[result]

(* {{1, 0.899155, 0.399155}, {1, 0.899155, 0.399155}} *)

For your example of a line surface intersection NMinimize cam get trapped in a local minimum with the default method:

par = 5 {Cos[u], Cos[v] + Sin[u], Sin[v]};
line = {1, 2, 3} + {4, -5, 6} t;
NMinimize[{SquaredEuclideanDistance[par, line], 0 < u < 2 Pi, -Pi < v < Pi}, {t, u, v}]
(* result: {2.97068, {t -> 0.39428, u -> 5.45283, v -> 0.964654}} *)

Notice the distance 2.97 is too high. Try a different method to get a better result:

NMinimize[{SquaredEuclideanDistance[par, line], 0 < u < 2 Pi, -Pi < v < Pi},
 {t, u, v}, Method -> "DifferentialEvolution"]
(* {1.77932*10^-16, {t -> -1.07503, u -> 2.29164, v -> -0.761533}} *)

Answered by flinty on November 27, 2020

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP