2
$\begingroup$

I am trying to use a system of equations to solve for the common tangent of a (polar) InterpolatingFunction and a circle. Unfortunately, both NSolve and FindRoot don't yield results, even when I "hand-feed" values into the equations into the equations (i.e. in combinations of equations 1,3,4 / 2,3,4 or 1,2,3)

My guess is that this may have to do with "excessive" precision, although I'm not sure.

Question: what function would I best use to get the closest approximations, if it is a precision problem? For my purposes, three decimal places would be sufficient.

(If someone sees a better basic geometric strategy to go at the problem, please do say.)

Please have a look at the graphic and the code:

diagram

(* intf[ang1] is the InterpolatingFunction (green),
rotated about {0,0} until it is tangent to the circle.
Goal: find out ang1, ang3, angrotation *)

intpoints = {{0, 17.9817}, {0.0521424, 18.5701}, {0.105454, 
    19.1892}, {0.15991, 19.8416}, {0.215494, 20.5299}, {0.272199, 
    21.2567}, {0.330029, 22.0243}, {0.388993, 22.8348}, {0.449111, 
    23.6897}, {0.510411, 24.5896}, {0.572931, 25.5338}, {0.636716, 
    26.5201}, {0.701823, 27.5438}, {0.76832, 28.5972}, {0.836286, 
    29.6684}, {0.905813, 30.7409}, {0.977011, 31.7919}, {1.05, 
    32.7921}, {1.12493, 33.7052}, {1.20195, 34.4887}, {1.28123, 
    35.0956}, {1.36295, 35.4775}, {1.44731, 35.5905}};
intf = Interpolation[intpoints]

(* values used in example *)
r1 = 9.5;
x1 = 34.2584;
y1 = 12;

(* 1. parametric equations for a circle of radius R1 with center \
located at {X1,Y1} *)
xcirclefn[ang2_] = r1*Cos[ang2] + x1;
ycirclefn[ang2_] = r1*Sin[ang2] + y1;

(* 2. equations for locating points based on rotation of axes of intf \
*)
xintf[ang1_] = intf[ang1]*Cos[ang1];
yintf[ang1_] = intf[ang1]*Sin[ang1];
xrot[ang1_, angrotation_] = 
  xintf[ang1]*Cos[angrotation] - yintf[ang1]*Sin[angrotation];
yrot[ang1_, angrotation_] = 
  xintf[ang1]*Sin[angrotation] + yintf[ang1]*Cos[angrotation];

(* 3. equations for angles of tangents *)
slopecircle[ang2_] = ycirclefn'[ang2]/xcirclefn'[ang2] ;
angleang3c[ang2_] = ArcTan[-1/slopecircle[ang2]] ;

slopeintf[ang1_] = yintf'[ang1]/xintf'[ang1];
angleang3i[ang1_, angrotation_] = 
  ArcTan[-1/slopeintf[ang1]] + angrotation;
(* angrotation is negative *)

(* Possible identities: *)
equation1 = xcirclefn[ang2] == xrot[ang1, angrotation];
equation2 = ycirclefn[ang2] == yrot[ang1, angrotation];
equation3 = angleang3c[ang2] == angleang3i[ang1, angrotation];
equation4 = intf[ang1] == Sqrt[xcirclefn[ang2]^2 + ycirclefn[ang2]^2];

(* CAD-measured values in radians:
ang1 = 0.714321 (CCW from intf[0.])
 ang2 = 3.09004 (CCW from 0.)
angrotation = -0.247318 (CW from 0.)
ang3 = -(2*Pi)+(Pi/2 + (-0.0515572)) = -4.76395 (CW from 0.) (1.51924 included angle) *)

(* Method 1:
NSolve[equation1&&equation2&&equation3,{ang1,ang2,angrotation}] *)

(* Method 2:
FindRoot[{equation1,equation3,equation4},{{ang1,0},{ang2,0},{\
angrotation,-Pi/2}}] *)
$\endgroup$
  • $\begingroup$ Could you complete your example with the interpolating function and the solving approaches you already tested? And maybe, if possible, switch to lowercase symbol names throughout? Uppercase sometimes collides with predefined symbols, causing all kinds of "surprising" behavior. $\endgroup$ – Jinxed Feb 22 '15 at 23:24
  • $\begingroup$ You might also want to define your functions using := (SetDelayed). $\endgroup$ – Jinxed Feb 22 '15 at 23:29
  • $\begingroup$ Depending on the Interpolation you used, there might not be a usable derivative. See the documentation for examples. $\endgroup$ – Jinxed Feb 22 '15 at 23:31
  • $\begingroup$ If you could provide me with your intf, I might be able to help. $\endgroup$ – Jinxed Feb 23 '15 at 0:42
  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 0) Browse the common pitfalls question 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – Dr. belisarius Feb 23 '15 at 6:58
2
$\begingroup$

Not sure if I understand, anyway:

intpoints = {{0, 17.9817}, {0.0521424, 18.5701}, {0.105454, 19.1892}, {0.15991, 19.8416}, 
            {0.215494, 20.5299}, {0.272199, 21.2567}, {0.330029, 22.0243}, {0.388993, 22.8348}, 
            {0.449111, 23.6897}, {0.510411, 24.5896}, {0.572931, 25.5338}, {0.636716, 26.5201}, 
            {0.701823, 27.5438}, {0.76832, 28.5972}, {0.836286,  29.6684}, {0.905813, 30.7409}, 
            {0.977011, 31.7919}, {1.05, 32.7921}, {1.12493, 33.7052}, {1.20195, 34.4887}, 
            {1.28123,  35.0956}, {1.36295, 35.4775}, {1.44731, 35.5905}};
intf = Interpolation[intpoints]

(*values used in example*)
r1 = 9.5; x1 = 34.2584; y1 = 12;
exp = RotationMatrix[-ang1].(intf[p] {Cos@p, Sin@p});
sol = NMinimize[{ang1, EuclideanDistance[{x1, y1}, exp] == r1 && 0 < p < Pi/2}, {p,  ang1}];
Show[ContourPlot[(x - x1)^2 + (y - y1)^2 == r1^2, {x, 0, 50}, {y, -10, 40}], 
     ParametricPlot[intf[p] {Cos@p, Sin@p}, {p, 0, Pi/2}], 
     ParametricPlot[exp /. sol[[2]], {p, 0.01, Pi/2}, PlotStyle -> Red]]

Mathematica graphics

| improve this answer | |
$\endgroup$
  • $\begingroup$ @Jinxed: info added, thanks! $\endgroup$ – jharrison Feb 23 '15 at 6:50
  • $\begingroup$ belisarius: thanks, will try adapting this method and see if I can pull the angles out. Looks much cleaner than mine. $\endgroup$ – jharrison Feb 23 '15 at 6:56
  • $\begingroup$ @jharrison There you have it done $\endgroup$ – Dr. belisarius Feb 23 '15 at 6:57
  • $\begingroup$ I started fresh and copied your code, and got: NMinimize::nsol: There are no points that satisfy the constraints {}. >> $\endgroup$ – jharrison Feb 23 '15 at 7:58
  • $\begingroup$ @jharrison: Same for me. $\endgroup$ – Jinxed Feb 23 '15 at 12:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.