1
$\begingroup$
ClearAll["Global`*"]
(*This program is to find the natural frquency of a beam with \
spring-mass connected at some arbitary location 'z'*)
Y = 2*10^11;
Iyy = 8.333*10^-6;
L = 4;
z = 0.01;
stiff = (Y*Iyy)/L^3;
k = {0.01, 1, 2, 4, 6, 8, 10, 25, 50, 100, 200, 300, 400, 500, 600, 
  800, 1*^3, 2*^3, 4*^3, 5*^3, 10*^3, 15*^3, 20*^3, 30*^3, 40*^3, 
  70*^3, 1*^5, 1*^6, 1*^8, 
  1*^12}(*THIS IS THE FIRST  PARAMETER REPRESENTS  STIFFNESS *)
m = {0.01, 1, 2, 4, 6, 8, 10, 25, 50, 100, 200, 300, 400, 500, 600, 
  800, 1*^3, 2*^3, 4*^3, 5*^3, 10*^3, 15*^3, 20*^3, 30*^3, 40*^3, 
  70*^3, 1*^5, 1*^6, 1*^8, 
  1*^12}(*THIS IS THE SECOND PARAMETER REPRESENTS  MASS*)
(*s[i_,j_]:=((m[[j]]*b^4)/(1-m[[j]]/k[[i]]*b^4))*)
(*Represents the scaling in right hand side of the differential \
equation *)
eq = (1/(2*b^3))*(((Sin[b*(1 - z)]*Sin[b*z])/
      Sin[b]) - ((Sinh[b*(1 - z)]*Sinh[b*z])/
      Sinh[b]));(*Green's function *)
eq1[i_, j_] := (1 - (((m[[j]]*b^4)/(1 - m[[j]]/k[[i]]*b^4))*eq));
sol = ArrayReshape[Table[eq1[i, j], {i, 30}, {j, 30}], {900, 1, 1}]
NSolve[sol[[1, 1]] == 0 && 0 < b < 10]

Now, this sol contains 900 equation how to solve for first 5 roots of this equation? How to solve all the equations contained in the sol in one single shot?

$\endgroup$
  • $\begingroup$ It's an eigenvalue problem, right? Then how about using NDEigensystem instead? $\endgroup$ – xzczd May 8 '18 at 8:07
  • $\begingroup$ Kind off, but when I used NSolve it works manually, but when I put in the loop it didn't . This is my try For[i = 1, i < 901, i++, S[[i]] = NSolve[sol[[i]] == 0 && 0 < b < 15]] $\endgroup$ – acoustics May 8 '18 at 8:11
  • $\begingroup$ Try FindInstance[...,5] instead of NSolve. I am busy to to go into deep of your code which is not completely correct. $\endgroup$ – user64494 May 8 '18 at 8:13
1
$\begingroup$

The equation solving in your loop is not stuck, it's just very slow in v11.2. You can speed up the code by removing the removable singularity in the equation:

Y = 2 10^11;
Iyy = (8333/1000)/10^6;
L = 4;
z = 1/100;
stiff = (Y Iyy)/L^3;
k = {1/100, 1, 2, 4, 6, 8, 10, 25, 50, 100, 200, 300, 400, 500, 600, 800, 1000, 2000, 
  4000, 5000, 10000, 15000, 20000, 30000, 40000, 70000, 100000, 1000000, 100000000, 
  1000000000000}
m = {1/100, 1, 2, 4, 6, 8, 10, 25, 50, 100, 200, 300, 400, 500, 600, 800, 1000, 2000, 
  4000, 5000, 10000, 15000, 20000, 30000, 40000, 70000, 100000, 1000000, 100000000, 
  1000000000000}
(* Notice the following modification. *)
eq = (Sin[b (1 - z)] Sin[b z] - Sin[b] (Sinh[b (1 - z)] Sinh[b z])/Sinh[b])/(2 b^3); 
eq1[i_, j_] := Sin[b] (1 - (m[[j]] b^4)/k[[i]]) - (m[[j]] b^4) eq;

test = Table[
    b /. NSolve[eq1[i, j] == 0 && 0 < b < 20, b, WorkingPrecision -> 32], {i, 30}, {j, 
     30}]; // AbsoluteTiming
(* {43.0112, Null} *)

Notice though the roots are quite close to the singularities, they're not identical, for example:

test[[1, 1, 1]] - 1
(* -8.220929635875291505012540*10^-8 *)

BTW, it's worth mentioning that v9 can solve these equations fast without the modification mentioned above, but it only finds 6252 roots, while v11.2 finds 6256.

$\endgroup$
  • $\begingroup$ ok, then from where the actual roots start then?, If it is throwing b->1 , there should be a zero cross over at 1 right? $\endgroup$ – acoustics May 8 '18 at 9:35
  • $\begingroup$ @vijay I have to retract my words, the roots are not spurious. Check my edit. $\endgroup$ – xzczd May 8 '18 at 12:44
  • $\begingroup$ ok, I will check, thanks for the help. $\endgroup$ – acoustics May 8 '18 at 16:23

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.