The acoustic problem of vibrations of bars is usually solved with graphic methods and most books that have these equations in them just mention that you could use the NewtonRaphson method to solve it. The approach in this article is different. This article uses the NuGet symbolic mathematics package MathNet.Symbolics together with my C# implementation of the JenkinTraub rootfinding algorithm of polynomials to solve these problems.
A Word of Warning
I might as well warn you from the start. If you don't know the acoustic equations that I'm solving here, this article is not going to be very useful for you. Or you might be a mathematician that has a similar math problem that requires solving an equation using a polynomial root finding method. For now, it's just a short article, and I might upgrade it later with some nice solution plots.
Introduction
The solution to the differential wave equation leads to a combination of two main functions with two different amplitudes, i.e., A*Sin
and B*Cos
. The boundary conditions imposed usually result in solving the equation for the (usually scaler) variables A and B. Typically, A and B are visualized as a standing incoming wave and leaving wave after D'Alembert's principle solution to the transverse wave equation. The problem is that the boundary conditions are often complex in nature, and there is no simple way of solving these equations, so generally, a numerical method is required. The question is just which one to use.
I will, for the most part, skip the derivations of how the equations are set up and deal with the solution straight away. All the solutions are solved using the continued fraction expansion and truncating it. This will lead to a Pade approximation to the function we are trying to solve. The denominator is removed by some simple algebraic manipulation such that a single complex polynomial rootfinding algorithm can be used. In this case, that is the JenkinsTraub algorithm. The implementation in pseudocode:

Replace function like:
\(\sin(x)\)
with a fraction:
\( \sin(x) = \frac{a_0 + a_1 \cdot x + a_2 \cdot x^2 + \cdots + a_n \cdot x^n}{b_0 + b_1 \cdot x + b_2 \cdot x^2 + \cdots + b_n \cdot x^n }\)

Rearrange it into a polynomial, i.e., remove all fractions and set the equation equal to zero
 Solve the truncated polynomial with a rootfinding technique.
The reason that the Pade approximation is so useful in acoustics is due to the solutions to the wave equation. They are often truncated alternating series, and it is in these cases that the Pade approximation or Shanks transform is most useful. Due to the alternating power series, the solutions tend either to infinity or negative infinity very fast due to the limitations of computational numbers.
For the transverse wave equations, a different approach is taken. The starting point (as always, by the way. Including the Pade approximation) is the Taylor series. Since the function here is two know Taylor series, they are just symbolically multiplied together, and the root algorithm is set to find the roots from this.
The physical difference between the longitudinal wave and the transverse wave can be illustrated quite simply. Assume that someone jumps at either one of the two locations:
The transfer of energy into the brown beam (bar) from jumps at the longitudinal section will result in compressions in the vertical direction. While jumping at the transverse location will result in a transverse wave along the grey beam (bar). The astute might realize that the words actually describe the difference between the direction of the initial force compared to the wave direction. The longitudinal wave moves in the same direction as the force is applied. The transverse wave travels at a 90degree angle of the applied force. While jumping up and down, the wave travels from side to side in the gray beam.
Yes I am simplifying a bit, since jumping in the longitudinal section will also result in a transverse wave, but the explanation is hopefully simpler to understand if I just mention this afterward. In solid structures, they almost always coincide, but not in fluids. Here longitudinal waves are the only possible ones (well, mostly. You can have special situations where they do play a primary role).
Longitudinal Wave Equation
Free Mass Loaded Bar
The equation to solve for is this:
$\begin{equation} \tan(k L) =  \frac{m}{m_b} k L \end{equation}$
The solution using Linq and the Symbolic package:
private void xCotxSolver(int N = 11)
{
List<string> ContFraction_xCotx = new List<string>();
ContFraction_xCotx.Add("1");
for (int i = 1; i <= N; i++)
{
ContFraction_xCotx.Add("x*x");
ContFraction_xCotx.Add((2 * (i + 1)  1).ToString());
}
double m_b = 1;
double m = 1;
ContinuedFraction xCotx = new ContinuedFraction(ContFraction_xCotx);
Expr equation = xCotx.Numerator + (m_b / m) * xCotx.Denomiator;
List<Complex> roots = RealPolynomialRoots.Roots
(equation.ReverseRealCoefficients("x"));
List<double> result =
roots.
FindAll(x => (x.Real >= 0)).
Select(x => x.Real).
ToList<double>();
result.Sort();
}
The solutions should be \( kL = 2.03, 4.91, 7.98, \cdots \).
The Freely Vibrating Bar  General Boundary Condition
The equation that one seeks the solution to:
$\begin{equation} \tan(k L) = j\frac{\left( Z_{m0} / \left(\rho_L c\right) \right) + \left( Z_m / \left( \rho_L c \right) \right)}{1+\left( Z_{m0} / \left(\rho_L c\right) \right)\left( Z_m / \left( \rho_L c \right) \right)} \end{equation}$
The code that does this:
private void FreeVibrationBar_GeneralBC()
{
double rho_L = 1;
double C = 1;
Complex Zm0 = new Complex(10, 1);
Complex Zm = new Complex(10, 1);
Complex u = new Complex(0, 1);
u *= (Zm0 / (rho_L * C) + (Zm / (rho_L * C))) /
(1 + (Zm0 / (rho_L * C)) * (Zm / (rho_L * C)));
List<string> ContFraction_Tan_x = new List<string>();
ContFraction_Tan_x.Add("0");
ContFraction_Tan_x.Add("x");
ContFraction_Tan_x.Add("1");
for (int i = 2; i < 25; i++)
{
ContFraction_Tan_x.Add("x*x");
ContFraction_Tan_x.Add((2 * (i)  1).ToString());
}
ContinuedFraction Tanx = new ContinuedFraction(ContFraction_Tan_x);
var TanXnumerator = Tanx.Numerator.Coefficients(Expr.Variable("x"));
var TanXDenominator = Tanx.Denomiator.Coefficients(Expr.Variable("x"));
var PolynomialAtanExpansion = (Tanx.Numerator  Tanx.Denomiator *
new MathNet.Numerics.Complex32((float)u.Real, (float)u.Imaginary)).
Coefficients(Expr.Variable("x")).
ToList().
Select(x => x.ComplexNumberValue);
List<Complex> roots_Atanx = ComplexPolynomialRoots.Roots
(PolynomialAtanExpansion.Reverse<Complex>().ToArray());
List<Complex> Positive_roots_Atanx =
roots_Atanx.
FindAll(xxr => (xxr.Real >= 0)).
Distinct<Complex>().ToList();
List<Complex> Atan = new List<Complex>();
int mm = 0;
foreach (var item in Positive_roots_Atanx)
{
Atan.Add(item);
Atan.Add(Math.PI * mm + Complex.Atan(u));
mm += 1;
}
}
Forced Vibrations on a Bar
The equation we seek the solution to is:
$\begin{equation} Z_{m0} = \rho_L c \frac{\left( Z_m / \rho_L c \right) + j \tan(k L)}{1 +\left( Z_m / \rho_L c \right) j \tan(k L) } \end{equation}$
By making the substitution, some algebra can be performed:
private void ForcedVibrationOnABar()
{
List<string> ContFraction_Tan_x = new List<string>();
ContFraction_Tan_x.Add("0");
ContFraction_Tan_x.Add("x");
ContFraction_Tan_x.Add("1");
for (int i = 2; i < 25; i++)
{
ContFraction_Tan_x.Add("x*x");
ContFraction_Tan_x.Add((2 * (i)  1).ToString());
}
ContinuedFraction Tanx = new ContinuedFraction(ContFraction_Tan_x);
double x1 = 10;
double r1 = 0;
var er1 = Expr.Parse(x1.ToString()).Multiply
(Tanx.Numerator.Multiply(Tanx.Numerator));
double val = r1 * r1 + x1 * x1  1;
var er2 = Expr.Parse((val).ToString()).Multiply
(Tanx.Numerator).Multiply(Tanx.Denomiator);
var er3 = Expr.Parse(x1.ToString()).Multiply(Tanx.Denomiator.Pow(2));
var tyty = er1
.Add(er2)
.Subtract(er3)
.RationalSimplify(Expr.Variable("x"));
var TanX = tyty
.Coefficients(Expr.Variable("x"))
.Select(x => x.RealNumberValue)
.ToList<double>()
.Reverse<double>()
.ToArray<double>();
List<Complex> roots_tanxx = RealPolynomialRoots.Roots(TanX);
var tw = roots_tanxx.
Select(x => x.Real).
ToList<double>().
FindAll(x => x > 0).
Distinct<double>().
ToList<double>();
tw.Sort();
List<double> tw2 = new List<double>();
for (int uuu = 0; uuu < tw.Count(); uuu++)
{
tw2.Add(Math.PI * uuu + 0.099669);
tw2.Add(Math.PI * uuu  1.4711);
}
tw2 = tw2.FindAll(x => x > 0).Distinct().ToList();
tw2.Sort();
List<double> test = new List<double>();
for (int hg = 0; hg < tw.Count; hg++)
{
test.Add(tw[hg]);
test.Add(tw2[hg]);
}
}
Transverse Wave Equation
The general solution to the transverse wave equation for a bar:
$\begin{equation} y(x,t) = \left[ A \cdot \cosh(g \cdot x) + B \cdot \sinh(g \cdot x) + C \cdot \cos(g \cdot x ) + D \cdot \sin(g \cdot x )\right] \cos (\omega t + \phi ) \end{equation}$
In order to solve this equation, we need its derivative, second derivative and third derivative:
$\begin{equation} \frac{d \, y(x,t)}{dx} = g \cdot \left[ A \cdot \sinh(g \cdot x) + B \cdot \cosh(g \cdot x)  C \cdot \sin(g \cdot x) + D \cdot \cos(g \cdot x) \right] \cos(\omega t + \phi ) \end{equation}$
Second derivative:
$\begin{equation} \frac{d^2 \, y(x,t)}{{dx}^2} = g^2 \cdot \left[A \cdot \cosh(g \cdot x) + B \cdot \sinh(g \cdot x)  C \cdot \cos(g \cdot x)  D \cdot \sin(g \cdot x) \right] \cos (\omega t + \phi ) \end{equation}$
Third derivative:
$\begin{equation} \frac{d^3 \, y(x,t)}{{dx}^3} =g^3 \cdot \left[A \cdot \sinh(g \cdot x) + B \cdot \cosh(g \cdot x) + C \cdot \sin(g \cdot x)  D \cdot \cos(g \cdot x) \right] \cos (\omega t + \phi ) \end{equation}$
These equations are used together with the appropriate boundary conditions:
Fixed end:
$\begin{equation} y = 0 \qquad \text{and} \qquad \frac{\partial y}{\partial x} = 0 \end{equation}$
Free end:
$\begin{equation} \frac{\partial^2 y}{{\partial x}^2} = 0 \qquad \text{and} \qquad \frac{\partial^3 y}{{\partial x}^3} = 0 \end{equation}$
Simply supported end:
$\begin{equation}y = 0 \qquad \text{and} \qquad \frac{\partial^2 y}{{\partial x}^2} = 0 \end{equation}$
Boundary conditions (BC) for the clamped end. This is usually the fixed beam since y = 0 and the derivative is also zero. This means that the initial position does not change at all, and the beam is always at a straight angle from its fixed point regardless of the vibrations present in the rest of the beam. The involves solving the equation:
$\begin{equation}\cos(g \cdot L) \cosh(g \cdot L) = 1 \end{equation}$
This is almost equal to the equation with a vibrating bar free at both ends. Think of this as the frequency that is generated by throwing a bar at the floor.
$\begin{equation}\cos(g \cdot L) \cosh(g \cdot L) = 1 \end{equation}$
With these two equations, the Pade approximation is not used (although you could). Instead, the Taylor (or rather the MacLaurin) series is used directly.
private void ClamedFreeBarVibration()
{
Expr xr = Expr.Parse("x");
Expr one = Expr.Parse("1");
Expr cosh = Expr.Parse("1");
Expr cos = Expr.Parse("1");
for (int t = 1; t <= 25; t++)
{
cosh += xr.Pow(2 * t) / Factorial(2 * t);
cos += Math.Pow(1, t) * xr.Pow(2 * t) / Factorial(2 * t);
}
Expr coshcos = cosh.Multiply(cos).Add(one).RationalSimplify(xr);
var Coefficients = coshcos
.Coefficients(Expr.Variable("x"))
.Select(x => x.RealNumberValue)
.ToList<double>()
.Reverse<double>()
.ToList<double>();
List<Complex> roots_coshcos = RealPolynomialRoots.Roots(Coefficients.ToArray());
List<Double> roots_cch =
roots_coshcos.
Where(p => ((p.Imaginary <= 0.01) && (p.Real >= 0))).
Select(x => x.Real).
Distinct<double>().
ToList<double>().
FindAll(x => Math.Abs(coshcos.Evaluate
(new Dictionary<string, FloatingPoint> { { "x", x } }).RealValue) < 1e4).
Select(x => Math.Round(x * 2 / Math.PI, 3)).
Distinct<double>().
ToList<double>();
roots_cch.Sort();
}
The simply supportedfree bar’s transcendental equation is:
$\begin{equation} \frac{\tan(g L)} {\tanh(gL)} = 1 \end{equation}$
Reference
Fundamentals of Acoustics, 4^{th} edition, Kinsler, Frey, Coppers and Sanders.
History
 20^{th} March, 2021: Initial release
 24^{th} March, 2021: Corrected written formula (code correct) for freely vibrating bar, added referances and fixed some minor formatting issues