How to find the roots of a very big equation. and why NSolve fails to extract the roots

ClearAll["Global*"];
Clear[b]
L = 1;
z[1] = L/4;
z[2] = (2*L)/4;
z[3] = (3*L)/4;
Y = 2*10^11;(*Youngs modulus *)
Iyy = 8.333*10^-6;(*area moment of inertia*)
A = 0.1^2;(*cross sectional area*)
mass = 1; (*mass at the free end of the cantielver*)
ρ = 7850;
n = 3; (* number of springs *)
W[1] = a[1]*Sin[b*x] + a[2]*Cos[b*x] + a[3]*Sinh[b*x] + a[4]*Cosh[b*x];
W[2] = a[5]*Sin[b*x] + a[6]*Cos[b*x] + a[7]*Sinh[b*x] + a[8]*Cosh[b*x];
W[3] = a[9]*Sin[b*x] + a[10]*Cos[b*x] + a[11]*Sinh[b*x] +
a[12]*Cosh[b*x];
W[4] = a[13]*Sin[b*(x - z[3])] + a[14]*Cos[b*(x - z[3])] +
a[15]*Sinh[b*(x - z[3])] + a[16]*Cosh[b*(x - z[3])];
w = Piecewise[{{W[1], x <= z[1]}, {W[2], z[1] <= x < z[2]}, {W[3],
z[2] <= x < z[3]}, {W[4], x > z[3]}}];

(*CANTILEVER BC*)
boundary[i_, j_] :=
Module[ {bc},
bc1 = {W[i] /. {x -> 0}, (D[W[i], {x, 1}]) /. {x -> 0}, (D[
W[j], {x, 2}]) /. {x -> L}, ((D[W[j], {x, 3}]) /. {x -> L})};
bc = bc1]

countinuity[i_, j_] :=
Module[{eq},
eq1 = {((W[i] /. x -> z[i]) - (W[j] /.
x -> z[i])), (((D[W[i], {x}]) /.
x -> z[i]) - ((D[W[j], {x}]) /.
x -> z[i])), (((D[W[i], {x, 2}]) /.
x -> z[i]) - ((D[W[j], {x, 2}]) /.
x -> z[i])), (((D[W[i], {x, 3}]) /.
x -> z[i]) - ((D[W[j], {x, 3}]) /. x -> z[i])) + (K[i]*
W[i] /. x -> z[i])}; eq = eq1 ]

e1 = boundary[1, 4];
e2 = countinuity[1, 2];
e3 = countinuity[2, 3];
e4 = countinuity[3, 4];
comb = Tuples[{0, 1*^12}, 3];
eq = Flatten[{e1, e2, e3, e4}];
var = Table[a[i], {i, 1, Length[eq]}];
R = Normal@CoefficientArrays[eq, var][[2]];
MatrixForm[R]
K[1] = 0; K[2] = 0; K[3] = 0;
P = Det[R]
s1 = P;
Plot[s1, {b, 0, 10}]
s2 = NSolve[s1 == 0 && 0 < b < 30]


I have a matrix whose dimension is 16 cross 16. I have take a determinant of that matix which leads to very huge equation and it is transcendental in nature. I wanted to find the roots of this equation, for that I used NSolve, but it is not giving me nay results. what are the other ways to extract the roots. and why NSolve fails for these kind of equation.

• Please do not include extraneous code. – Michael E2 Nov 22 '18 at 2:39
• Change the definition of s1 to s1 = P // Simplify; Note the difference in the complexity of the expressions; LeafCount /@ {P, s1} evaluates to {7590, 12}. In addition to improving efficiency it will also improve precision of calculations. Compare the precision for P /. b -> 10.15 // FullForm with s1 /. b -> 10.15 // FullForm – Bob Hanlon Nov 22 '18 at 4:29

If you restrict the range of b NSolve can solve your problem

NSolve[{s1 == 0, 1 < b < 10}, b]
{{b -> 1.8751}, {b -> 4.69409}, {b -> 7.85476}}

s1 = P // Simplify

(* 128 b^24 (1 + Cos[b] Cosh[b]) *)


The exact solutions are expressed as Root objects

(s2 = Solve[{s1 == 0, 0 < b < 30}]) // Column


The approximate numerical roots are

s2 // N

(* {{b -> 1.8751}, {b -> 4.69409}, {b -> 7.85476}, {b -> 10.9955}, {b ->
14.1372}, {b -> 17.2788}, {b -> 20.4204}, {b -> 23.5619}, {b ->
26.7035}, {b -> 29.8451}} *)

• Bob, This is even faster. I never tried using Solve for a transcendental equation in Mathematica. Thanks a lot. Actually, the thing is K[1], K[2], K[3] are the three defined variables. these variables take two values 0 and infinity. If I take a combination of these variables there are 8 such combinations. I have written a code but seems like it running forever. can you give me some input to implement this properly. – acoustics Nov 22 '18 at 6:08
• Nice. Bursts my balloon, though. This was the first time I recall NSolve finding multiple roots at a simple root because of numerical noise. But now it seems unnecessary. – Michael E2 Nov 22 '18 at 13:15
• If any K[n] is Infinity then Det[R] will be infinite and your problem cannot be solved as posed. Perhaps you mean K[n] can be extremely large and you want to look at the behavior in that case. Let the non-zero value be the variable inf. Then the term with the highest order of 'inf' will dominate. Perhaps what you want is NSolve[{CoefficientList[Det[R /. kValues] // FullSimplify, inf][[-1]] == 0, 0 < b < 30}, b]. This will in general be very slow. – Bob Hanlon Nov 22 '18 at 15:12

As @BobHanlon has shown there is a "simple" solution to the specific problem. As for the OP's question "Why," I was taken in by the numerical behavior of s1, the inability of WorkingPrecision in NSolve to resolve the problem, and the trick I used to post-process the NSolve results.

The trouble doesn't seem to have to do with the numerics of evaluating s1 per se, although 16-28 digits are lost for b ranging from 14.1 to 23.6 and the derivative at the roots range from 10^35 to 10^45. This is where the trouble begins and is likely the principal cause of the failure of NSolve; however, evaluating s1 at high precision produces smoothly changing values, and FindRoot has no trouble finding the roots. The trouble must lie in the internal algorithm used by NSolve, which is unknown (to me).

The "trick" is to effectively use a lower "PrecisionGoal" by reducing the precision of the result returned by NSolve and removing duplicate roots. NSolve does not have PrecisionGoal and AccuracyGoal options. It does have an option Method -> {"Tolerance" -> <number(s)>}, but it's not clear that it can be applied here. I'll leave the answer, since I do not think this approach has been shown on this site.

Higher working precision helps, especially for b > 8 (it is obvious in the plot that round-off error dominates at greater values of b). In this case, numerical noise in the last few digits, at whatever WorkingPrecsion, means several roots show up multiple times at slightly different locations. Lowering the final precision makes these roots the same.

s2 = DeleteDuplicates@ N@ NSolve[s1 == 0 && 0 < b < 30, WorkingPrecision -> 32]
(*
{{b -> 1.8751},  {b -> 4.69409}, {b -> 7.85476}, {b -> 10.9955},
{b -> 14.1372}, {b -> 17.2788}, {b -> 20.4204}, {b -> 23.5619}}
*)

• In a way, if I define the working precision I can get the roots. Is that what you are saying – acoustics Nov 22 '18 at 4:32
• Changing the definition of s1 to s1 = P // Simplify enables use of machine precision. s2 = NSolve[s1 == 0 && 0 < b < 30] and will give the ten roots in the interval {0, 30} – Bob Hanlon Nov 22 '18 at 5:14
• But there are some roots which get repeated, and I used DeleteDuplicates[b/.s2], but could not able to remove the duplicates. – acoustics Nov 22 '18 at 5:18
• With v11.3 on my Mac and the revised definition of s1 (i.e., including Simplify), s2 = NSolve[s1 == 0 && 0 < b < 30] evaluates to {{b -> 1.8751}, {b -> 4.69409}, {b -> 7.85476}, {b -> 10.9955}, {b -> 14.1372}, {b -> 17.2788}, {b -> 20.4204}, {b -> 23.5619}, {b -> 26.7035}, {b -> 29.8451}}` There are no duplicate roots. – Bob Hanlon Nov 22 '18 at 5:21