# Unable to find root for user-defined function

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:

However, when I attempt to solve

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


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.

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]


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]


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

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

four possible solutions cy for given cx.

• 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). Nov 12, 2019 at 15:55
• Thanks for that code Ulrich. I'll look into explaining it. Nov 12, 2019 at 17:29

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}]


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

• The definition of getPointLen2[t1, t2] is missing! Nov 13, 2019 at 7:18
• Ok sorry. I had renamed it to getSolution. Updated code above to reflect change. Nov 13, 2019 at 11:57

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"]


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;
(* 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"}]


• Thanks Tim. That's great! Now that's a new kind of geometry the ancients would envy :) I need to upgrade to 12. Nov 16, 2019 at 12:18
• 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. Nov 16, 2019 at 14:03