# Unable to solve an equation using NSolve

I have a Euler Bernoulli beam supported by a linear spring at L/3, and mass at L/2, and a torsional spring at 0.75L. I am interested in finding the natural frequencies of this system. Using general beam solution I constructed the piecewise filed variable w which gives the deformed shape for a particular frequency. after applying boundary condition and compatibility condition across springs and mass I have formed the matrix R, The det of R will gives the frequencies. The det equation is P, I am interested in finding the roots of this equation. I used NSolve to get those roots. If Find the null space associated with the first root, using the NUllspace function of matrix R, and If I plot the equation w from 0 to L, The plot should look something like this. But I am not getting this. ClearAll["Global*"];
Clear[b]

L = 4;
z1 = L/3;
z2 = L/2;
z3 = (3*L)/4;

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

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

e = (w1 /. {x -> z1}) - (w2 /. {x -> z1});
e = ((D[w1, {x}]) /. {x -> z1}) - ((D[w2, {x}]) /. {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});

e = (w2 /. {x -> z2}) - (w3 /. {x -> z2});
e = ((D[w2, {x}]) /. {x -> z2}) - ((D[w3, {x}]) /. {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 -> z1});

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][];
P = Det[R];

KT = 10^12; M = 10^12; KR = 10^12;
s1 = P;
s2 = NSolve[s1 == 0 && 0 < b < 10]
s3 = N[b /. s2]
s4 = s3[];
{uu, ww, vv} = SingularValueDecomposition[R /. b -> s4];
NN = Last[Transpose[vv]];   sub1 = Flatten[{var, b}]; sub2 =
Flatten[{NN, s4}];    m =
w /. Table[sub1[[i]] -> sub2[[i]], {i, 1, Length[sub1]}]
Plot[m, {x, 0, L}]

• As I understand the code works, but you did not get what you wanted? – Alex Trounev Feb 10 at 15:31

There is an error in the boundary conditions that I fixed. But other options are possible. I used FindRoot[], which is optimal in this case.

ClearAll["Global*"];
Clear[b]

L = 4;
z1 = L/3;
z2 = L/2;
z3 = (3*L)/4;

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

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

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});

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*(w2 /. {x -> z1});

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][];
P = Det[R];

KT = 10^12; M = 10^12; KR = 10^12;
s1 = P;

s2 = Table[
FindRoot[s1 == 0, {b, n}, WorkingPrecision -> 100,
MaxIterations -> 1000], {n, 1.5, 10, 1}]

s3 = N[b /. s2]
Table[s4 = s3[[i]];
{uu, ww, vv} = SingularValueDecomposition[R /. b -> s4];
NN = Last[Transpose[vv]]; sub1 = Flatten[{var, b}];
sub2 = Flatten[{NN, s4}];
m = w /. Table[sub1[[i]] -> sub2[[i]], {i, 1, Length[sub1]}];
Plot[m, {x, 0, 4}, PlotRange -> Automatic,
PlotLabel -> Row[{"b=", s3[[i]]}]], {i, 1, Length[s3]}] • Trounev, thanks solved my problem, can you please tell what was the problem with boundary conditions – acoustics Feb 11 at 5:18
• There is a small glitch the equation e, the last term M*(w2 /. {x -> z1}) is supposed to be this M*(w2 /. {x -> z2}). there is a typo from my side in the OP. This creates the mess for the first root plot. – acoustics Feb 11 at 6:14
• @acoustics Check e, e in my and in your code. – Alex Trounev Feb 11 at 12:43
• Actually, I copied the code from my two spring case, I didn't look at the boundary conditions. my fault. Thanks Alex – acoustics Feb 11 at 15:41
• @acoustics you're welcome! – Alex Trounev Feb 11 at 16:41

If you restrict the range of b NSolve can solve the equation (after some time)

NSolve[{s1 == 0 && 1 < b < 10}, b, Reals]
(*{{b -> 1.43604}, {b -> 2.35619}, {b -> 2.9357}, {b -> 3.96197},
{b -> 4.71239}, {b -> 7.06858}, {b -> 7.24871}}*)

• I used GraphicsMeshFindIntersections function to find the roots, and later I found the null space, which clearly violates by boundary conditions. I also got the same set of roots using the graphics mesh function. I don't know what mistake I am making which violates the physics of the problem . – acoustics Feb 10 at 11:03
• You didn't show your version of MeshFindIntersections , without knowledge of this attempt it's hard to find the mistake! – Ulrich Neumann Feb 10 at 15:27
• NSolve will give more accurate results if you simplify s1 first: s1 = P // Simplify // ComplexExpand[#, TargetFunctions -> {Re, Im}] & // Simplify; – Bob Hanlon Feb 10 at 17:43
• s1 = P // FullSimplify; is better yet. – Bob Hanlon Feb 10 at 18:32

You are looking for eigenvalues? It should be widespread knowledge that finding the roots of the characteristic polynomial is the worst numerical way to compute eigenvalues. Why don't you use tools that where specifically designed for this task, for example Eigenvalues?

• Ok, could you please brief me how to set up eigenvalue problem for the system I described in OP – acoustics Feb 10 at 11:04
• Ah, now I see that R depends nonlinearly on b, so my suggestion is not that easy to implement. Nevertheless: Have you plotted Abs[P] over the range of possible b? This function is quite a beast! What I try to say: The current algorithm for computing the frequencies is not adequate. – Henrik Schumacher Feb 10 at 11:16
• You should consider reformulating your problem into an eigenvalue problem. Usually, "frequencies" of interest are related to eigenfunction of a suitable linear operator. In this case, this operator could be related to the Hessian of the energy, evaluated at the configuration of interest. But your question is actually a bit to vague to be sure. – Henrik Schumacher Feb 10 at 11:16
• I did not understand the later energy part of the configuration, what does it mean? – acoustics Feb 10 at 13:01
• Typically the notion of "natural frequencies" of equilibrium state of a dynamical system correspond to eigenfunctions of the Hessian of the total energy of the system. More precisely, the eigenfunctions correspond to generalized(!) eigenvectors with respect to the inner product induced by mass. I cannot be more concrete, because I would have to set up the whole mechnical system, discretize it, find the static equilibria, compute the Hessian and so on. – Henrik Schumacher Feb 10 at 13:49