1
$\begingroup$

I created the function getSolution[xVal,yVal] for a geometry problem I'm working on and wish to solve it for getSolution[2,x]==0. I can plot it clearly showing the roots: enter image description here

However, when I attempt to solve

 FindRoot[getSolution[2, x] == 0, {x, 1.5}]

I receive the errors:

Part::partd: Part specification x$18835[[2]] is longer than depth of object.

Part::partd: Part specification x$18835[[1]] is longer than depth of object.

Which I've tracked down I believe to FindRoot attempting to obtain the derivative of the function for the Newton method. Since if I code:

D[getSolution[2, x], x]

I receive the same message. I was wondering if someone could offer some suggestions to fix this code? Here's my code:

getAngle[a_, b_, c_] := Module[{line1, line2, line3, aVal, bVal, val},
   aVal = a - b;
   bVal = c - b;
   val = (aVal.bVal)/(EuclideanDistance[{0, 0}, 
        aVal] EuclideanDistance[{0, 0}, bVal]);
   ArcCos[val]
   ];

makeEqn[p1_, p2_] := 
  Function[{x}, (p2[[2]] - p1[[2]])/(p2[[1]] - p1[[1]]) (x - 
       p1[[1]]) + p1[[2]]];

getLength[a_, b_] := EuclideanDistance[a, b];

getSolution[cx_, cy_] := 
 Module[{c, m, acLength, macAngle, d, ebLength, bmLength, emLength, 
   acEqn, acAngle, e, deEqn, cbEqn, f, x, a, b},
  a = {0, 0};
  b = {7, 0};
  c = {cx, cy};
  m = {getLength[a, b]/2, 0};
  acLength = getLength[a, c];
  macAngle = getAngle[m, a, c];
  d = {acLength Cos[macAngle], 0};

  bmLength = getLength[b, m];
  acEqn = makeEqn[a, c];

  e = {x, acEqn[x]} /. 
     NSolve[acEqn[x]^2 + (x - m[[1]])^2 == bmLength^2 && x > 0, x] // 
    First;
  deEqn = makeEqn[d, e];
  cbEqn = makeEqn[c, b];
  f = {x, deEqn[x]} /. NSolve[deEqn[x] == cbEqn[x], x] // First;
  getLength[e, b] - getLength[f, b] // N
  ]

Not sure this is needed but here's the problem: Essentially need to find the point C in the diagram such that the triangle FEB is isosceles. When getSolution[x,y]=0, it is. Also, don't wish to just generate a table of values and find it manually for one case as above as I wish to solve the larger problem getSolution[x,y]=0 for all x and y in applicable domain.

enter image description here

$\endgroup$
2
$\begingroup$

In fact, for at least your specific case, NSolve[] is not really necessary, as one can explicitly solve the underlying linear and quadratic equations:

getSolution[cx_, cy_] := Module[{a, b, c, m, acLength, d, bmLength, e, f, x},
   a = {0, 0}; b = {7, 0}; c = {cx, cy}; m = {EuclideanDistance[a, b]/2, 0};
   acLength = EuclideanDistance[a, c];
   d = {acLength Normalize[m - a].Normalize[c - a], 0};
   bmLength = EuclideanDistance[b, m];
   e = (c (cx m[[1]] + Sqrt[(bmLength Norm[c])^2 - cy^2 m[[1]]^2]))/Norm[c]^2;
   f = LinearSolve[{Cross[e - d], Cross[b - c]}, {(d - e).Cross[d], (c - b).Cross[c]}];
   EuclideanDistance[e, b] - EuclideanDistance[f, b]]

ContourPlot[Evaluate[getSolution[cx, cy] == 0 /. Abs -> Identity],
            {cx, 0, 3}, {cy, -4, 4}, FrameLabel -> Automatic, PlotPoints -> 75]

solution curve


One can in fact do some further analysis to separate out the two branches of the solution. Skipping the details,

ContourPlot[{-7 cx + cx^2 + cy^2 == 0,
             -343 cx + 147 cx^2 - 21 cx^3 + cx^4 + 147 cy^2 - 21 cx cy^2 - cy^4 == 0},
            {cx, 0, 3}, {cy, -4, 4}, FrameLabel -> Automatic, PlotPoints -> 75]

solution branches

$\endgroup$
  • $\begingroup$ NSolve works with your code! Try NSolve[{getSolution[1, cy] == 0, 0 < cy}, cy] (*{{cy -> 1.31843}, {cy -> 2.44949}, {cy -> 11.1473}}*) $\endgroup$ – Ulrich Neumann Nov 13 at 14:06
1
$\begingroup$

Not an answer, but perhaps a small step further on:

You try to solve the equation getSolution[cx, cy]==0 .

First I changed the definition of getSolution[cx, cy] (thanks to @John Doty 's comment)

getSolution[cx_?NumericQ, cy_?NumericQ] := 
Block[{c, m, acLength, macAngle, d, ebLength, bmLength, emLength, acEqn, acAngle, e, deEqn, cbEqn, f, x, a, b}, a = {0, 0};b = {7, 0};c = {cx, cy};m = {getLength[a, b]/2, 0};acLength = getLength[a, c];macAngle = getAngle[m, a, c];d = {acLength Cos[macAngle], 0};bmLength = getLength[b, m];acEqn = makeEqn[a, c];e = {x, acEqn[x]} /. NSolve[acEqn[x]^2 + (x - m[[1]])^2 ==bmLength^2 && x > 0, x] //First;deEqn = makeEqn[d, e];cbEqn = makeEqn[c,b];f = {x, deEqn[x]} /. NSolve[deEqn[x] == cbEqn[x], x] // First;
getLength[e, b] - getLength[f, b]]

If I now try to apply Contourplot to visualize the solution

ContourPlot[getSolution[cx, cy] == 0 , {cx, 0, 3}, {cy, -4, 4}, FrameLabel -> Automatic, Evaluated -> True]

Mathematica v12 gives an error message Part::partd: Part specification x$2040162[[2]] is longer than depth of object. which indicates that's something wrong with the function definition.

Surprisingly

ContourPlot[getSolution[cx, cy] , {cx, 0, 3}, {cy, -4, 4},Contours -> {0}, FrameLabel -> Automatic]

evaluates and shows

enter image description here

four possible solutions cy for given cx.

$\endgroup$
  • 2
    $\begingroup$ I believe the partd error arises from an attempt to symbolically analyze a procedure. This can be avoided by giving two starting values for x rather than one, or by restricting getSolution to numeric inputs by defining it as getSolution[cx_?NumericQ, cy_?NumericQ]:=... (be sure to clear its old definition first). $\endgroup$ – John Doty Nov 12 at 15:55
  • $\begingroup$ Thanks for that code Ulrich. I'll look into explaining it. $\endgroup$ – Dominic Nov 12 at 17:29
1
$\begingroup$

Thanks guys. The ?NumericQ did it. Here's the plot in 3D with the yellow line representing the points (x,y) which produce the isosceles triangle:

mycontour[t1_] := t2 /. FindRoot[getSolution[t1, t2], {t2, 1.4}];
pp1 = ParametricPlot3D[{t1, mycontour[t1], 0}, {t1, 1, 3}, 
   PlotStyle -> Yellow];
pp2 = Plot3D[getSolution[t, t2], {t, 1, 3}, {t2, 1, 3}, 
   BoxRatios -> {1, 1, 1}, AxesLabel -> {"x", "y", "z"}];
Show[{pp2, pp1}]

enter image description here

And here's a (numerically) exact solution in a Manipulate environement: enter image description here

$\endgroup$
  • $\begingroup$ The definition of getPointLen2[t1, t2] is missing! $\endgroup$ – Ulrich Neumann Nov 13 at 7:18
  • $\begingroup$ Ok sorry. I had renamed it to getSolution. Updated code above to reflect change. $\endgroup$ – Dominic Nov 13 at 11:57
1
$\begingroup$

If you have triangles inscribed in circles, there is a good chance that synthetic geometry and FindGeometricConjectures will find some simplifying assumptions given your constraints.

The following code suggest that if FEB is isosceles, then triangle CEB is a right triangle.

gsc[cx_, cy_] := Module[{aa, bb, ri},
  aa = {0, 0};
  bb = {7, 0};
  ri = RandomInstance[
    GeometricScene[
     {a -> aa, b -> bb, c -> {cx, cy}, d, e, f, 
      r}, {p1 == Style[Polygon[{b, f, e}], Cyan], 
      EuclideanDistance[b, f] == EuclideanDistance[b, e], 
      EuclideanDistance[r, f] < EuclideanDistance[r, c], 
      GeometricAssertion[{c, f, r, b}, "Collinear"], 
      GeometricAssertion[{d, f, e}, "Collinear"], 
      GeometricAssertion[{a, c, e}, "Collinear"], 
      CircleThrough[{b, e, c, d}, r]}
     ], RandomSeeding -> 4
    ];
  {ri, FindGeometricConjectures[ri]}
  ]
{ri, fc} = gsc[1.75, 1.5]
ri["Points"]

Conjectures

Although it did not fall out from the conjectures, another useful relation is the Intersecting Chords Theorem that can be verified using the points:

{Norm[c - f] Norm[f - b], Norm[d - f] Norm[f - e]} /. 
 ri["Points"]
(* {4.12068, 4.12068} *)

Another constraint to note is that the solution is only valid for $\measuredangle ACB\geq \frac{\pi }{2}$, which you can use Solve to find the relations:

Last[
 cy /. Simplify[
   Solve[{π/2 == 
      Last@PolygonAngle@Triangle[{{0, 0}, {bx, 0}, {cx, cy}}]}, 
    cy], {bx, cx} ∈ Reals]]
First[cx /. 
  Simplify[Solve[{π/2 == 
      Last@PolygonAngle@Triangle[{{0, 0}, {bx, 0}, {cx, cy}}]}, 
    cx], {bx, cy} ∈ Reals]]
(* Sqrt[(bx - cx) cx] *)
(* 1/2 (bx - Sqrt[bx^2 - 4 cy^2]) *)

Knowing the point C defines the center and radius of the circle. You can use RegionIntersection to solve for E. All other points can be determined from E. Below is a simple demonstration.

ep[cx_, cy_, bx_: 7] := 
 Module[{cf, sal, green, pf, a, b, c, d, r, radius, circ, rc, rl, 
   rint, o, e, ebc, f, tri},
  cf = RGBColor[0.39215686274509803`, 0.5843137254901961, 
    0.9294117647058824];
  sal = RGBColor[
   0.9803921568627451, 0.5019607843137255, 0.44705882352941173`];
  green = RGBColor[
   0.19607843137254902`, 0.803921568627451, 0.19607843137254902`];
  pf = {#3, Disk[#1, 0.25], White, Text[#2, #1, Background -> #3]} &;
  a = {0, 0};
  c = {cx, cy};
  b = {bx, 0};
  (* Circle Center and Radius *)
  r = (b + c)/2;
  radius = Norm[b - c]/2;
  circ = Circle[r, radius];
  (* Find point E by Circle Line Intersection *)
  rc = Region[circ];
  rl = InfiniteLine[a, c];
  rint = RegionIntersection[rc, rl][[1]][[1]];
  o = Ordering[Norm[#] & /@ rint];
  e = rint[[Last@o]];
  (* Rotate BE to align with CB to create F *)
  ebc = Pi/2 - PlanarAngle[c -> {b, e}];
  f = RotationTransform[ebc, b][e];
  tri = Triangle[{e, f, b}];
  (* Use Intersecting Chord Theorem to  create D *)
  d = f + Norm[c - f] Norm[f - b] Normalize[f - e]/Norm[f - e];
  Graphics[{Line[{a, b}], Opacity[0.25], 
    EdgeForm[Directive[Thick, Blue]], Opacity[1], Thick, Line[{c, a}],
     Line[{c, e}], Line[{c, f}], Line[{d, f}], Cyan, tri, Black, 
    EdgeForm[Directive[Thin, Black]], circ, pf[a, "A", sal], 
    pf[b, "B", sal], pf[c, "C", green], pf[e, "E", cf], 
    pf[r, "R", cf], pf[f, "F", cf], pf[d, "D", cf]}, 
   PlotRange -> {{-0.5, bx + 0.8}, {-3 bx/7, 2/3 bx}}]
  ]
Manipulate[
 If[cy > 0.99 Sqrt[(bx - cx) cx], cy = 0.99 Sqrt[(bx - cx) cx]]; 
 ep[cx, cy, bx], {{cx, 2}, 1/2 (bx - Sqrt[bx^2 - 4 cy^2]), 4, 
  Appearance -> "Labeled"}, {{cy, 1.5}, 0.001, 
  0.99 Sqrt[(bx - cx) cx], Appearance -> "Labeled"}, {{bx, 7}, 5, 10, 
  Appearance -> "Labeled"}]

Manipulate Demo

$\endgroup$
  • $\begingroup$ Thanks Tim. That's great! Now that's a new kind of geometry the ancients would envy :) I need to upgrade to 12. $\endgroup$ – Dominic Nov 16 at 12:18
  • 1
    $\begingroup$ Thanks Dominic. I updated the post with a constrained Manipulate that I think should run on earlier versions. I am intrigued by the new functionality. It has a CAD and GeoGebra vibe to it that can help one identify constraints that may not be obvious. $\endgroup$ – Tim Laska Nov 16 at 14:03

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.