Unable to plot the determinant function (Possible bug in CoefficientArrays?)

Question

ClearAll["Global`*"]
(*Data*)
L1 = 4;
L2 = 4;
Iyy1 = (0.1*0.1^4)/12;
Iyy2 = (0.1*0.1^4)/12;
A1 = 0.1*0.1;
A2 = 0.1*0.1;
ρ1 = 7850;
ρ2 = 7850;
Y1 = 2*10^11;
Y2 = 2*10^11;
b1 = Surd[(ρ1*A1*ω^2*L1^4)/(Y1*Iyy1), 4];
b2 = Surd[(ρ1*A2*ω^2*L2^4)/(Y2*Iyy2), 4];
λ1 = Sqrt[(A1*L1^2)/Iyy1];
λ2 = Sqrt[(A2*L2^2)/Iyy2];
v1 = FullSimplify[Surd[b1^4/λ1^2, 4]];
v2 = FullSimplify[Surd[b2^4/λ2^2, 4]];
(*Beam Functions*)
W1 = FullSimplify[
   C1*Cos[b1*x1] + C2*Sin[b1*x1] + C3*Cosh[b1*x1] + C4*Sinh[b1*x1]];
W2 = FullSimplify[
   C5*Cos[b2*x2] + C6*Sin[b2*x2] + C7*Cosh[b2*x2] + C8*Sinh[b2*x2]];
(*Bar function*)
U1 = FullSimplify[C9*Cos[v1*x1] + C10*Sin[v1*x1]];
U2 = FullSimplify[C11*Cos[v2*x2] + C12*Sin[v2*x2]];

(*Boundary condition*)
e1 = FullSimplify[W1 /. x1 -> 0];
e2 = FullSimplify[(D[W1, {x1, 1}]) /. x1 -> 0];
e3 = FullSimplify[U1 /. x1 -> 0];
e4 = FullSimplify[W2 /. x2 -> 0];
e5 = FullSimplify[(D[W2, {x2, 1}]) /. x2 -> 0];
e6 = FullSimplify[U2 /. x2 -> 0];

(*Compatability condition*)
(*Displacement contunity*)
e7 = FullSimplify[(W1 /. x1 -> L1) - (U2 /. x2 -> L2)];
e8 = FullSimplify[(W2 /. x2 -> L2) + (U1) /. x1 -> L1];
(*Slope Contunity*)
e9 = FullSimplify[((D[W1, {x1, 1}]) /. 
      x1 -> L1) - ((D[W2, {x2, 1}]) /. x2 -> L2)];
(*Moment Contunity*)
e10 = FullSimplify[((Y1*Iyy1)/L1*((D[W1, {x1, 2}]) /. x1 -> L1)) + ((
      Y2*Iyy2)/L2*((D[W2, {x2, 2}]) /. x2 -> L2))];
(*Force Contuinity*)
e11 = FullSimplify[((Y1*Iyy1)/
      L1^2*((D[W1, {x1, 3}]) /. x1 -> L1)) - (Y2*
      A2 ((D[U2, {x2, 1}]) /. x2 -> L2))];
e12 = FullSimplify[((Y2*Iyy2)/
      L2^2*((D[W2, {x2, 3}]) /. x2 -> L2)) + (Y1*
      A1 ((D[U1, {x1, 1}]) /. x1 -> L1))];
(*Solving*)
R = FullSimplify[
   Normal@CoefficientArrays[{e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, 
       e11, e12}, {C1, C2, C3, C4, C5, C6, C7, C8, C9, C10, C11, 
       C12}][[2]]];
R1 = MatrixForm[R];
MatrixRank[R];
P = FullSimplify[Det[R]]
Plot[P, {ω, 0, 500}]
s1 = NSolve[P == 0 && 0 < ω < 500]
s2 = Flatten[ω /. s1];
s3 = s2[[i]];
fn = s3/(2*π)

I have a 12 homogenous equation. I have written these equations in MatrixForm. And I have taken the Det of that matrix R which happened to be the function of Omega, I tried to plot that function but I am getting the following error Integer expected at position 2. So tried the Following 1. Used 'FullSimplify' 2.Expand 3.TrigReduce 4.TrigExpand But still, I could not able to solve the issue. I guess maybe I have used Surd which leads to the fourth root in the determinant function. How to overcome this error. This error is stopping me from solving for the roots of determinant function

Answer 1

Minimal example that reproduces error

Fails

Normal@CoefficientArrays[{Times[1., y, Surd[x, 2]]}, y]
(* {{0.}, {{1. Surd[x, 2.]}}} *)

These do not fail

Normal@CoefficientArrays[Times[1., y, Surd[x, 2]], y]
(* {0, {1. Surd[x, 2]}} *)

Normal@CoefficientArrays[{Times[y, Surd[x, 2]]}, y]
(* {{0}, {{Surd[x, 2]}}} *)

It may be a bug, but for whatever reason CoefficientArrays is transforming

Surd[x, 2]

into

Surd[x, 2.]

And the evaluation fails because the second arguments is expected to be an integer, and obviously

IntegerQ[2.]
(* False *)

Solutions

After

R = FullSimplify[
   Normal@CoefficientArrays[{e1, e2, e3, e4, e5, e6, e7, e8, e9, e10, 
       e11, e12}, {C1, C2, C3, C4, C5, C6, C7, C8, C9, C10, C11, 
       C12}][[2]]];

you can patch the result

R = (R /. Surd[a_, b_] :> Surd[a, Round[b]])

Or use Block as sugested by @CarlWoll in the comments.

Block[
 {Surd},
 Normal@CoefficientArrays[
   {Times[1., y, Surd[x, 2]]}
   , y]
 ]

Or, see below the suggestion by Wolfram Support

ClearAttributes[Surd,NumericFunction]

Reply from Wolfram Support

This is the first reply I got for [CASE:4086017]

Thank you for contacting Wolfram Technical Support. I understand that you found a problem in the way Mathematica handles the arguments in Normal@CoefficientArrays. In general, the introduction of finite precision numbers converts everything to finite precision and only if infinite precision is given everywhere than the result is given with infinite precision. The problem that you encountered is generated by the presence of the float number 1. and all coefficients are floating points If integer 1 would be used instead the problem would not appear

Normal@CoefficientArrays[{Times[1,y,Surd[x,2]]},y]

returns

{{0},{{Surd[x, 2]}}}

Use of Rationalize function would render all these coefficients as integers

Rationalize[Normal@CoefficientArrays[{Times[1.,y,Surd[x,2]]},y]]

gives

{{0.},{{1. Surd[x, 2.]}}}

Best regards,

Wolfram Technical Support Wolfram Research Inc.

And then on further communication

Thank you for contacting Wolfram Technical Support. Thank you for clarifying the issue. Indeed the presence of a finite precision arguments transforms all arguments into finite precision ones leading to unintended consequence of changing all arguments to lowest precision. This includes arguments that must be integers. Surd[] has NumericFunction as Attribute. Clearing it prevents this type of argument type changing.

Normal@CoefficientArrays[{Times[1.,y,Surd[x,2]]},y]/.{x->0.5}

returns an error

Surd::int: Integer expected at position 2 in Surd[0.5, 2.]. {{0.}, {{1. Surd[0.5, 2.]}}}

while

ClearAttributes[Surd,NumericFunction]
Normal@CoefficientArrays[{Times[1.,y,Surd[x,2]]},y]/.{x->0.5}

returns the numerical evaluations

{{0.},{{0.707107}}}

I filed a report to our developers to consider this.