This model describes string vibration where tension increases with elongation:
\[\rho \text{A} \frac{\partial^{2}}{\partial t^{2}} y{\left(x,t \right)} = \left(\text{T}_0 + \frac{\text{A} \text{E}}{2 \text{L}} \int\limits_{0}^{\text{L}} \left(\frac{\partial}{\partial x} y{\left(x,t \right)}\right)^{2}\, dx\right) \frac{\partial^{2}}{\partial x^{2}} y{\left(x,t \right)} \]
Where:
Replacing derivatives with finite differences:
\[\rho \text{A} \left(- \frac{2 y{\left(x,t \right)}}{{\Delta}t^{2}} + \frac{y{\left(x,t - {\Delta}t \right)}}{{\Delta}t^{2}} + \frac{y{\left(x,t + {\Delta}t \right)}}{{\Delta}t^{2}}\right) = \left(\frac{\text{A} \text{E} \int\limits_{0}^{\text{L}} \left(- \frac{y{\left(x - {\Delta}x,t \right)}}{2 {\Delta}x} + \frac{y{\left(x + {\Delta}x,t \right)}}{2 {\Delta}x}\right)^{2}\, dx}{2 \text{L}} + \text{T}_0\right) \left(- \frac{2 y{\left(x,t \right)}}{{\Delta}x^{2}} + \frac{y{\left(x - {\Delta}x,t \right)}}{{\Delta}x^{2}} + \frac{y{\left(x + {\Delta}x,t \right)}}{{\Delta}x^{2}}\right) \]
Applying discrete indices (\(t \rightarrow \text{n} {\Delta}t\), \(x \rightarrow \text{i} {\Delta}x\)) and approximating the integral using a summation over the grid points:
\[\frac{\rho \text{A} \left({y}_{\text{i},\text{n} + 1} + {y}_{\text{i},\text{n} - 1} - 2 {y}_{\text{i},\text{n}}\right)}{{\Delta}t^{2}} = \left(\frac{\text{A} \text{E} {\Delta}x \sum_{\text{j}=0}^{\text{M} - 1} \left(\frac{{y}_{\text{j} + 1,\text{n}}}{2 {\Delta}x} - \frac{{y}_{\text{j} - 1,\text{n}}}{2 {\Delta}x}\right)^{2}}{2 \text{L}} + \text{T}_0\right) \left(\frac{{y}_{\text{i} + 1,\text{n}}}{{\Delta}x^{2}} + \frac{{y}_{\text{i} - 1,\text{n}}}{{\Delta}x^{2}} - \frac{2 {y}_{\text{i},\text{n}}}{{\Delta}x^{2}}\right) \]
Solving for the next time step yields the update scheme:
\[{y}_{\text{i},\text{n} + 1} = - {y}_{\text{i},\text{n} - 1} + 2 {y}_{\text{i},\text{n}} + \frac{{\Delta}t^{2} \left(\frac{\text{A} \text{E} {\Delta}x \sum_{\text{j}=0}^{\text{M} - 1} \left(\frac{{y}_{\text{j} + 1,\text{n}}}{2 {\Delta}x} - \frac{{y}_{\text{j} - 1,\text{n}}}{2 {\Delta}x}\right)^{2}}{2 \text{L}} + \text{T}_0\right) \left(\frac{{y}_{\text{i} + 1,\text{n}}}{{\Delta}x^{2}} + \frac{{y}_{\text{i} - 1,\text{n}}}{{\Delta}x^{2}} - \frac{2 {y}_{\text{i},\text{n}}}{{\Delta}x^{2}}\right)}{\rho \text{A}} \]
Or, organizing the update scheme:
\[{y}_{\text{i},\text{n} + 1} = \left(\text{T}_0 + {\text{Z}}_{2} \sum_{\text{j}=0}^{\text{M} - 1} \left({y}_{\text{j} + 1,\text{n}} - {y}_{\text{j} - 1,\text{n}}\right)^{2}\right) \left({y}_{\text{i} + 1,\text{n}} + {y}_{\text{i} - 1,\text{n}} - 2 {y}_{\text{i},\text{n}}\right) {\text{Z}}_{1} - {y}_{\text{i},\text{n} - 1} + 2 {y}_{\text{i},\text{n}} \]
Where:
\[{\text{Z}}_{1} = \frac{{\Delta}t^{2}}{\rho \text{A} {\Delta}x^{2}} \]
\[{\text{Z}}_{2} = \frac{\text{A} \text{E}}{8 \text{L} {\Delta}x} \]
Analyzing the stability of this nonlinear scheme is complex. We can approximate it by linearizing the nonlinear tension term. Assume the dominant contribution to the integral comes from a characteristic displacement shape.
We approximate the squared slope term within the integral. Consider a simple case like a triangular pluck shape with maximum amplitude \(\alpha \text{L}\) at the midpoint \(\frac{\text{L}}{2}\). The slope magnitude would be approximately \(2 \alpha = 2\alpha\).
We replace the term \(\left(\frac{{y}_{\text{j} + 1,\text{n}}}{2 {\Delta}x} - \frac{{y}_{\text{j} - 1,\text{n}}}{2 {\Delta}x}\right)^{2}\) inside the summation with a constant approximation based on a characteristic amplitude ratio \(\alpha\):
\[\left(\frac{{y}_{\text{j} + 1,\text{n}}}{2 {\Delta}x} - \frac{{y}_{\text{j} - 1,\text{n}}}{2 {\Delta}x}\right)^{2} \approx \left(\frac{\alpha \text{L}}{{\Delta}x}\right)^{2} = \frac{\alpha^{2} \text{L}^{2}}{{\Delta}x^{2}} \]
Here, \(\alpha\) represents a typical ratio of maximum amplitude to string length.
Substituting this approximation into the update scheme:
\[{y}_{\text{i},\text{n} + 1} = - {y}_{\text{i},\text{n} - 1} + 2 {y}_{\text{i},\text{n}} + \frac{{\Delta}t^{2} \left(\frac{\alpha^{2} \text{A} \text{E} \text{L} \text{M}}{2 {\Delta}x} + \text{T}_0\right) \left(\frac{{y}_{\text{i} + 1,\text{n}}}{{\Delta}x^{2}} + \frac{{y}_{\text{i} - 1,\text{n}}}{{\Delta}x^{2}} - \frac{2 {y}_{\text{i},\text{n}}}{{\Delta}x^{2}}\right)}{\rho \text{A}} \]
Now we substitute the wave-form ansatz \({y}_{\text{i},\text{n}} = \text{Y} e^{i \left(\text{i} k {\Delta}x + \text{n} w {\Delta}t\right)}\) into the linearized update scheme:
\[\text{Y} e^{i \left(\text{i} k {\Delta}x + w {\Delta}t \left(\text{n} + 1\right)\right)} = 2 \text{Y} e^{i \left(\text{i} k {\Delta}x + \text{n} w {\Delta}t\right)} - \text{Y} e^{i \left(\text{i} k {\Delta}x + w {\Delta}t \left(\text{n} - 1\right)\right)} + \frac{{\Delta}t^{2} \left(\frac{\alpha^{2} \text{A} \text{E} \text{L} \text{M}}{2 {\Delta}x} + \text{T}_0\right) \left(- \frac{2 \text{Y} e^{i \left(\text{i} k {\Delta}x + \text{n} w {\Delta}t\right)}}{{\Delta}x^{2}} + \frac{\text{Y} e^{i \left(\text{n} w {\Delta}t + k {\Delta}x \left(\text{i} - 1\right)\right)}}{{\Delta}x^{2}} + \frac{\text{Y} e^{i \left(\text{n} w {\Delta}t + k {\Delta}x \left(\text{i} + 1\right)\right)}}{{\Delta}x^{2}}\right)}{\rho \text{A}} \]
Expanding terms:
\[\text{Y} e^{i w {\Delta}t} e^{i \text{i} k {\Delta}x} e^{i \text{n} w {\Delta}t} = \frac{\alpha^{2} \text{E} \text{L} \text{M} \text{Y} {\Delta}t^{2} e^{i k {\Delta}x} e^{i \text{i} k {\Delta}x} e^{i \text{n} w {\Delta}t}}{2 \rho {\Delta}x^{3}} - \frac{\alpha^{2} \text{E} \text{L} \text{M} \text{Y} {\Delta}t^{2} e^{i \text{i} k {\Delta}x} e^{i \text{n} w {\Delta}t}}{\rho {\Delta}x^{3}} + \frac{\alpha^{2} \text{E} \text{L} \text{M} \text{Y} {\Delta}t^{2} e^{- i k {\Delta}x} e^{i \text{i} k {\Delta}x} e^{i \text{n} w {\Delta}t}}{2 \rho {\Delta}x^{3}} + 2 \text{Y} e^{i \text{i} k {\Delta}x} e^{i \text{n} w {\Delta}t} - \text{Y} e^{- i w {\Delta}t} e^{i \text{i} k {\Delta}x} e^{i \text{n} w {\Delta}t} + \frac{\text{T}_0 \text{Y} {\Delta}t^{2} e^{i k {\Delta}x} e^{i \text{i} k {\Delta}x} e^{i \text{n} w {\Delta}t}}{\rho \text{A} {\Delta}x^{2}} - \frac{2 \text{T}_0 \text{Y} {\Delta}t^{2} e^{i \text{i} k {\Delta}x} e^{i \text{n} w {\Delta}t}}{\rho \text{A} {\Delta}x^{2}} + \frac{\text{T}_0 \text{Y} {\Delta}t^{2} e^{- i k {\Delta}x} e^{i \text{i} k {\Delta}x} e^{i \text{n} w {\Delta}t}}{\rho \text{A} {\Delta}x^{2}} \]
Dividing by the common factor \(\text{Y} e^{i \text{i} k {\Delta}x} e^{i \text{n} w {\Delta}t}\) and rearranging:
\[e^{i w {\Delta}t} + e^{- i w {\Delta}t} = \frac{\alpha^{2} \text{E} \text{L} \text{M} {\Delta}t^{2} e^{i k {\Delta}x}}{2 \rho {\Delta}x^{3}} - \frac{\alpha^{2} \text{E} \text{L} \text{M} {\Delta}t^{2}}{\rho {\Delta}x^{3}} + \frac{\alpha^{2} \text{E} \text{L} \text{M} {\Delta}t^{2} e^{- i k {\Delta}x}}{2 \rho {\Delta}x^{3}} + 2 + \frac{\text{T}_0 {\Delta}t^{2} e^{i k {\Delta}x}}{\rho \text{A} {\Delta}x^{2}} - \frac{2 \text{T}_0 {\Delta}t^{2}}{\rho \text{A} {\Delta}x^{2}} + \frac{\text{T}_0 {\Delta}t^{2} e^{- i k {\Delta}x}}{\rho \text{A} {\Delta}x^{2}} \]
Simplifying using Euler’s formula to get the dispersion relation for the linearized system:
\[2 \cos{\left(w {\Delta}t \right)} = \frac{\alpha^{2} \text{A} \text{E} \text{L} \text{M} {\Delta}t^{2} \left(\cos{\left(k {\Delta}x \right)} - 1\right) + 2 \rho \text{A} {\Delta}x^{3} + 2 \text{T}_0 {\Delta}t^{2} {\Delta}x \left(\cos{\left(k {\Delta}x \right)} - 1\right)}{\rho \text{A} {\Delta}x^{3}} \]
Solving for the angular frequency \(w\):
\[w = \left[ \frac{- \operatorname{acos}{\left(\frac{\alpha^{2} \text{A} \text{E} \text{L} \text{M} {\Delta}t^{2} \cos{\left(k {\Delta}x \right)} - \alpha^{2} \text{A} \text{E} \text{L} \text{M} {\Delta}t^{2} + 2 \rho \text{A} {\Delta}x^{3} + 2 \text{T}_0 {\Delta}t^{2} {\Delta}x \cos{\left(k {\Delta}x \right)} - 2 \text{T}_0 {\Delta}t^{2} {\Delta}x}{2 \rho \text{A} {\Delta}x^{3}} \right)} + 2 \pi}{{\Delta}t}, \ \frac{\operatorname{acos}{\left(\frac{\alpha^{2} \text{A} \text{E} \text{L} \text{M} {\Delta}t^{2} \cos{\left(k {\Delta}x \right)} - \alpha^{2} \text{A} \text{E} \text{L} \text{M} {\Delta}t^{2} + 2 \rho \text{A} {\Delta}x^{3} + 2 \text{T}_0 {\Delta}t^{2} {\Delta}x \cos{\left(k {\Delta}x \right)} - 2 \text{T}_0 {\Delta}t^{2} {\Delta}x}{2 \rho \text{A} {\Delta}x^{3}} \right)}}{{\Delta}t}\right] \]
For \(w\) to be real, the argument of arccosine must be between -1 and 1. We also substitute the number of spatial steps \(\text{M} \approx \frac{\text{L}}{{\Delta}x}\) (approximating floor):
\[-1 \le \frac{\alpha^{2} \text{A} \text{E} \text{L}^{2} {\Delta}t^{2} \cos{\left(k {\Delta}x \right)} - \alpha^{2} \text{A} \text{E} \text{L}^{2} {\Delta}t^{2} + 2 \rho \text{A} {\Delta}x^{4} + 2 \text{T}_0 {\Delta}t^{2} {\Delta}x^{2} \cos{\left(k {\Delta}x \right)} - 2 \text{T}_0 {\Delta}t^{2} {\Delta}x^{2}}{2 \rho \text{A} {\Delta}x^{4}} \le 1 \]
\[- 4 \rho \text{A} {\Delta}x^{4} \le - \alpha^{2} \text{A} \text{E} \text{L}^{2} {\Delta}t^{2} - 2 \text{T}_0 {\Delta}t^{2} {\Delta}x^{2} + \left(\alpha^{2} \text{A} \text{E} \text{L}^{2} {\Delta}t^{2} + 2 \text{T}_0 {\Delta}t^{2} {\Delta}x^{2}\right) \cos{\left(k {\Delta}x \right)} \le 0 \]
The right inequality \(- \alpha^{2} \text{A} \text{E} \text{L}^{2} {\Delta}t^{2} - 2 \text{T}_0 {\Delta}t^{2} {\Delta}x^{2} + \left(\alpha^{2} \text{A} \text{E} \text{L}^{2} {\Delta}t^{2} + 2 \text{T}_0 {\Delta}t^{2} {\Delta}x^{2}\right) \cos{\left(k {\Delta}x \right)} \le 0\) is always met. The left inequality \(- 4 \rho \text{A} {\Delta}x^{4} \le - \alpha^{2} \text{A} \text{E} \text{L}^{2} {\Delta}t^{2} - 2 \text{T}_0 {\Delta}t^{2} {\Delta}x^{2} + \left(\alpha^{2} \text{A} \text{E} \text{L}^{2} {\Delta}t^{2} + 2 \text{T}_0 {\Delta}t^{2} {\Delta}x^{2}\right) \cos{\left(k {\Delta}x \right)}\) provides the stability condition.
Substituting the worst-case value \(\cos{\left(k {\Delta}x \right)} = -1\) gives the most restrictive condition:
\(- 4 \rho \text{A} {\Delta}x^{4} \le - 2 \alpha^{2} \text{A} \text{E} \text{L}^{2} {\Delta}t^{2} - 4 \text{T}_0 {\Delta}t^{2} {\Delta}x^{2}\)
Solving for \({\Delta}x\) gives the stability condition for the linearized Kirchhoff-Carrier scheme:
\({\Delta}x \ge \frac{\sqrt{2} \sqrt{\text{T}_0 {\Delta}t^{2} + {\Delta}t \sqrt{2 \alpha^{2} \rho \text{A}^{2} \text{E} \text{L}^{2} + \text{T}_0^{2} {\Delta}t^{2}}}}{2 \sqrt{\rho} \sqrt{\text{A}}}\)
Substituting the expression for tension based on the fundamental frequency:
\[\frac{\sqrt{2} \sqrt{\text{L}} \sqrt{{\Delta}t} \sqrt{\sqrt{2} \sqrt{\rho} \sqrt{\alpha^{2} \text{E} + 8 \rho \text{L}^{2} f_{0}^{4} {\Delta}t^{2}} + 4 \rho \text{L} f_{0}^{2} {\Delta}t}}{2 \sqrt{\rho}} \le {\Delta}x \]
\[\sqrt{\frac{\text{L} {\Delta}t \left({\Delta}t f_{0}^{2} \cdot 4 \rho \text{L} + \sqrt{2 \rho \left(\alpha^{2} \text{E} + {\Delta}t^{2} f_{0}^{4} \cdot 8 \rho \text{L}^{2}\right)}\right)}{2 \rho}} \leq {\Delta}x \]
\[\frac{\sqrt{\text{L}} \frac{\sqrt{2}}{2 \sqrt{\rho}} \sqrt{\sqrt{2} \sqrt{\rho} \sqrt{\alpha^{2} \text{E} + \frac{8 \rho \text{L}^{2} f_{0}^{4}}{f_{s}^{2}}} + \frac{4 \rho \text{L} f_{0}^{2}}{f_{s}}}}{\sqrt{f_{s}}} \leq \frac{\text{L}}{\text{M}} \]
\[\frac{4 \rho \text{L} f_{0}}{\text{M}} \leq \sqrt{- 2 \alpha^{2} \rho \text{E} + \frac{4 \rho^{2} \text{L}^{2} f_{s}^{2}}{\text{M}^{4}}} \]
\[f_{0} \leq \frac{\text{M} \sqrt{- 2 \alpha^{2} \rho \text{E} + \frac{4 \rho^{2} \text{L}^{2} f_{s}^{2}}{\text{M}^{4}}}}{4 \rho \text{L}} \]
When \(f_{s} = 44100\), \(\text{M} = 3\), \(\text{L} = 1\), \(\rho = 7850\), \(\alpha = 0.01\), \(\text{E} = 200000000000.0\):
\[f_{0} \leq 7349.80501496208 \]
