# I have a working code, but work very slowly. want to improve the code

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


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).

• Did you plot the final solution? is it coming. Feb 4, 2019 at 13:42
• I am still struggling to get the desired results Feb 5, 2019 at 16:16