0
$\begingroup$

I have working code which gives the result as expected by me. But the problem is it is taking too much of time to give me those results. I have a determinant function which depends on z1,z2,z3, kt,m, KR,b. Out of these variables first six variables I know, the only unknown is b. I will input z1,z2,z3, kt,m, KR and find the roots of determinant which is a function of b alone now. For this, I have written two modules f and g. f takes input kt, m, kr and beta(beta root reader, if beta=1, means it as to read the first root.). After this, I have defined a variable called comb which consists of all possible combination of kt, m, kr (these three variables can take either, 0 or infinity). The other module is g, what it does is it takes this combination as input and throw my final result. after this I am just plotting my final result.

ClearAll["Global`*"];
Clear[b]
L = 1;
z1 = L/3;
z2 = L/2;
z3 = (3*L)/4;
n = 3;
w1 = a[1]*Sin[b*x] + a[2]*Cos[b*x] + a[3]*Sinh[b*x] + a[4]*Cosh[b*x];
w2 = a[5]*Sin[b*(x - z1)] + a[6]*Cos[b*(x - z1)] + 
   a[7]*Sinh[b*(x - z1)] + a[8]*Cosh[b*(x - z1)];
w3 = a[9]*Sin[b*(x - z2)] + a[10]*Cos[b*(x - z2)] + 
   a[11]*Sinh[b*(x - z2)] + a[12]*Cosh[b*(x - z2)];
w4 = a[13]*Sin[b*(x - z3)] + a[14]*Cos[b*(x - z3)] + 
   a[15]*Sinh[b*(x - z3)] + a[16]*Cosh[b*(x - z3)];

w = Piecewise[{{w1, x <= z1}, {w2, z1 <= x <= z2}, {w3, 
     z2 <= x <= z3}, {w4, x >= z3}}];

(*SS BC*)
e[1] = w1 /. {x -> 0};
e[2] = (D[w1, {x, 2}]) /. {x -> 0};
e[3] = w4 /. {x -> L};
e[4] = D[w4, {x, 2}] /. {x -> L};

(*Compatability condition for translation spring*)
e[5] = (w1 /. {x -> z1}) - (w2 /. {x -> z1});
e[6] = ((D[w1, {x, 1}]) /. {x -> z1}) - ((D[w2, {x, 1}]) /. {x -> z1});
e[7] = ((D[w1, {x, 2}]) /. {x -> z1}) - ((D[w2, {x, 2}]) /. {x -> z1});
e[8] = ((D[w1, {x, 3}]) /. {x -> z1}) - ((D[w2, {x, 3}]) /. {x -> 
       z1}) + kt*(w1 /. {x -> z1});

(*Forming matrix for translational springs *)
e[9] =  (w2 /. {x -> z2}) - (w3 /. {x -> z2});
e[10] = ((D[w2, {x, 1}]) /. {x -> z2}) - ((D[w3, {x, 1}]) /. {x -> 
       z2});
e[11] = ((D[w2, {x, 2}]) /. {x -> z2}) - ((D[w3, {x, 2}]) /. {x -> 
       z2});
e[12] = ((D[w2, {x, 3}]) /. {x -> z2}) - ((D[w3, {x, 3}]) /. {x -> 
       z2}) + m*(w1 /. {x -> z2});


(*Forming matrix for translational springs *)
e[13] = (w3 /. {x -> z3}) - (w4 /. {x -> z3});
e[14] = ((D[w3, {x, 1}]) /. {x -> z3}) - ((D[w4, {x, 1}]) /. {x -> 
       z3});
e[15] = ((D[w3, {x, 2}]) /. {x -> z3}) - ((D[w4, {x, 2}]) /. {x -> 
       z3}) + kr*((D[w3, {x, 1}]) /. {x -> z3});
e[16] = ((D[w3, {x, 3}]) /. {x -> z3}) - ((D[w4, {x, 3}]) /. {x -> 
       z3});

eq = Table[e[i], {i, 1, 16}];
var = Table[a[i], {i, 1, 16}];
R = Normal@CoefficientArrays[eq, var][[2]];
MatrixForm[R];
P = Expand[Det[R]];

f[KT_, M_, KR_, beta_] := 
  Module[{mm}, kt = KT; m = M; kr = KR; r = beta; s1 = P; 
   s2 = NSolve[{s1 == 0, 0 < b < 20}, b]; s3 = N[b /. s2]; 
   s4 = s3[[r]]; {uu, ww, vv} = 
    SingularValueDecomposition[R /. b -> s4]; 
   NN = Last[Transpose[vv]];   sub1 = Flatten[{var, b}]; 
   sub2 = Flatten[{NN, s4}];    
   mm = w /. Table[sub1[[i]] -> sub2[[i]], {i, 1, Length[sub1]}]; 
   Return[mm]];
comb = Tuples[{0, 1*^12}, 3]
g[i_, r_] := 
 Module[{s5}, spring = comb[[i]]; a1 = spring[[1]]; a2 = spring[[2]]; 
  a3 = spring[[3]]; a4 = r; s5 = f[a1, a2, a3, a4]]
modes = Table[g[i, 1], {i, 1, 2^n}];
Table[Plot[modes[[i]], {x, 0, L}, PlotRange -> All], {i, 1, 
  Length[modes]}]
$\endgroup$
2
$\begingroup$

As a partial solution I'd start with a case distinction instead of comb, and taking actual limits to infinity in order to simplify the expressions:

G[0, 0, 0] = FullSimplify[P /. {kt -> 0, m -> 0, kr -> 0}];
G[0, 0, ∞] = Limit[P/kr /. {kt -> 0, m -> 0}, kr -> ∞] // FullSimplify;
G[0, ∞, 0] = Limit[P/m /. {kt -> 0, kr -> 0}, m -> ∞] // FullSimplify;
G[0, ∞, ∞] = Limit[P/(m kr) /. {kt -> 0}, {m -> ∞, kr -> ∞}] // FullSimplify;
G[∞, 0, 0] = Limit[P/kt /. {m -> 0, kr -> 0}, kt -> ∞] // FullSimplify;
G[∞, 0, ∞] = Limit[P/(kt kr) /. {m -> 0}, {kt -> ∞, kr -> ∞}] // FullSimplify;
G[∞, ∞, 0] = Limit[P/(kt m) /. {kr -> 0}, {kt -> ∞, m -> ∞}] // FullSimplify;
G[∞, ∞, ∞] = Limit[P/(kt m kr), {kt -> ∞, m -> ∞, kr -> ∞}] // FullSimplify;

These are now easier to check for roots. If you have a look at these G[kt,m,kr] expressions you'll see that they factorize somewhat, and you can look for zeros of the factors.

For example, let's look at G[0, 0, ∞]. A scaled plot like

Plot[E^(-b)/b^21*G[0, 0, ∞], {b, 0, 30}]

shows zeros around $b=4,6,10,14,17,19,23,26,30$. Using these as starting points for a root search:

FindRoot[E^(-b)/b^21*G[0, 0, ∞] == 0, {b, #}] & /@ {4,6,10,14,17,19,23,26,30}
(* {{b -> 4.21587}, {b -> 6.28319}, {b -> 10.0041},
    {b -> 13.9317}, {b -> 16.902},  {b -> 18.8496},
    {b -> 22.5725}, {b -> 26.4987}, {b -> 29.4686}} *)

(there's also $b=0$ as a solution, of course).

$\endgroup$
  • $\begingroup$ Did you plot the final solution? is it coming. $\endgroup$ – acoustics Feb 4 at 13:42
  • $\begingroup$ I am still struggling to get the desired results $\endgroup$ – acoustics Feb 5 at 16:16

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.