# 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*Sin[b*x] + a*Cos[b*x] + a*Sinh[b*x] + a*Cosh[b*x];
w2 = a*Sin[b*(x - z1)] + a*Cos[b*(x - z1)] +
a*Sinh[b*(x - z1)] + a*Cosh[b*(x - z1)];
w3 = a*Sin[b*(x - z2)] + a*Cos[b*(x - z2)] +
a*Sinh[b*(x - z2)] + a*Cosh[b*(x - z2)];
w4 = a*Sin[b*(x - z3)] + a*Cos[b*(x - z3)] +
a*Sinh[b*(x - z3)] + a*Cosh[b*(x - z3)];

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

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

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

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

(*Forming matrix for translational springs *)
e = (w3 /. {x -> z3}) - (w4 /. {x -> z3});
e = ((D[w3, {x, 1}]) /. {x -> z3}) - ((D[w4, {x, 1}]) /. {x ->
z3});
e = ((D[w3, {x, 2}]) /. {x -> z3}) - ((D[w4, {x, 2}]) /. {x ->
z3}) + kr*((D[w3, {x, 1}]) /. {x -> z3});
e = ((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][];
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[]; a2 = spring[];
a3 = spring[]; 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. – acoustics Feb 4 '19 at 13:42
• I am still struggling to get the desired results – acoustics Feb 5 '19 at 16:16