| Hosted by CoCalc | Download
Kernel: Python 2 (SageMath)

Постановка задачи

Дан прямой однородный стержень. Для него заданы граничные условия:

  • dY(x)dxx=0=d3Y(x)dx3x=0=0\left . \frac{dY(x)}{dx} \right \vert_{x = 0} = \left . \frac{d^3Y(x)}{dx^3} \right \vert_{x = 0} =0

  • d2Y(x)dx2x=l=d3Y(x)dx3x=l=0\left . \frac{d^2Y(x)}{dx^2} \right \vert_{x = l} = \left . \frac{d^3Y(x)}{dx^3} \right \vert_{x = l} =0

Найти собственные значения и собственные формы деформации для поперечных колебаний прямого однородного стержня при отсутствии сжимающей нагрузки и трения в системе при заданных краевых условиях. Найти первые три нормированные собственные формы и построить их графики.

Поиск нормализованных собственных форм

Уравнение, описывающее изгибные колебания стержня, выглядит слеудющим образом:

EId4y(x,t)dx4+pd2y(x,t)dx2+εdy(x,t)dx+md2y(x,t)dt2=0EI \frac{d^4y(x, t)}{dx^4} + p \frac{d^2y(x, t)}{dx^2} + \varepsilon \frac{dy(x, t)}{dx} + m \frac{d^2y(x, t)}{dt^2} = 0

В силу отсутствия сжимающей нагрузки ε\varepsilon и вязкого трения pp уравнение изгибных колебаний принимает следующий вид:

EId4y(x,t)dx4+md2y(x,t)dt2=0EI \frac{d^4y(x, t)}{dx^4} + m \frac{d^2y(x, t)}{dt^2} = 0

Применим метод разделения переменных и найдём решения уравнения в виде y(x,t)=Y(x)T(t)y(x, t) = Y(x) \cdot T(t):

EIYIV(x)T(t)+mY(x)T(t)=0EI \cdot Y^{IV}(x) T(t) + m \cdot Y(x) T''(t) = 0EImYIV(x)Y(x)=T(t)T=λ\frac{EI}{m} \cdot \frac{Y^{IV}(x)}{Y(x)} = -\frac{T''(t)}{T} = \lambdaEIYIV(x)=λmy(x)EI Y^{IV}(x) = \lambda m y(x)YIV(x)λmEIy(x)=0Y^{IV}(x) - \frac{\lambda m}{EI} y(x) = 0

Сделаем замену параметра:β4=λmEI,  β>0\beta^4 = \frac{\lambda m}{EI}, \; \beta > 0.

Решаем характеристическое уравнение:

σ4β4=0\sigma^4 - \beta^4 = 0σk=βe2πk4i,    k=0,1,2,3\sigma_{k} = \beta e^{\frac{2\pi k}{4}i}, \;\; k = 0,1,2,3

Общее решение этого линейного однородного уравнения имеет вид:

Y(x)=C1sinβx+C2cosβx+C3sinhβx+C4coshβxY(x) = C_1 \sin \beta x + C_2 \cos \beta x + C_3 \sinh \beta x + C_4 \cosh \beta x

Вычислим некоторые производные общего решения в концах отрезка:

Y(x)=β(C1cosβxC2sinβx+C3coshβx+C4sinhβx)Y'(x) = \beta ( C_1 \cos \beta x - C_2 \sin \beta x + C_3 \cosh \beta x + C_4 \sinh \beta x )Y(0)=β(C1+C3)Y'(0) = \beta ( C_1 + C_3 )Y(x)=β2(C1sinβxC2cosβx+C3sinhβx+C4coshβx)Y''(x) = \beta^2 ( -C_1 \sin \beta x - C_2 \cos \beta x + C_3 \sinh \beta x + C_4 \cosh \beta x )Y(l)=β2(C1sinβlC2cosβl+C3sinhβl+C4coshβl)Y''(l) = \beta^2 ( -C_1 \sin \beta l - C_2 \cos \beta l + C_3 \sinh \beta l + C_4 \cosh \beta l )Y(x)=β3(C1cosβx+C2sinβx+C3coshβx+C4sinhβx)Y'''(x) = \beta^3 ( -C_1 \cos \beta x + C_2 \sin \beta x + C_3 \cosh \beta x + C_4 \sinh \beta x )Y(0)=β3(C1+C3)Y'''(0) = \beta^3 ( -C_1 + C_3 )Y(l)=β3(C1cos(βl)+C2sin(βl)+C3cosh(βl)+C4sinh(βl))Y'''(l) = \beta^{3} \bigl(- C_{1} \cos{\left (\beta l \right )} + C_{2} \sin{\left (\beta l \right )} + C_{3} \cosh{\left (\beta l \right )} + C_{4} \sinh{\left (\beta l \right )}\bigr)

С учётом граничных условий получается следующая система уравнений:

C1+C3=0,C_1 + C_3 = 0,C1sinβlC2cosβl+C3sinhβl+C4coshβl=0,-C_1 \sin \beta l - C_2 \cos \beta l + C_3 \sinh \beta l + C_4 \cosh \beta l = 0,C1+C3=0,-C_1 + C_3 = 0,C1cosβl+C2sinβl+C3coshβl+C4sinhβl=0-C_1 \cos \beta l + C_2 \sin \beta l + C_3 \cosh \beta l + C_4 \sinh \beta l = 0

Из этой системы уравнений сразу получается, что C1=C3=0C_1 = C_3 = 0, и в общем-то мы имеем дело с системой

C2cosβl+C4coshβl=0,-C_2 \cos \beta l + C_4 \cosh \beta l = 0,C2sinβl+C4sinhβl=0C_2 \sin \beta l + C_4 \sinh \beta l = 0

Если мы ищем нетривиальные решения, то её определитель должен равняться нулю:

cosβlcoshβlsinβlsinhβl=cosβlsinhβlcoshβlsinβl\left \vert \begin{array}{cc} -\cos \beta l & \cosh \beta l \\ \sin \beta l & \sinh \beta l \end{array} \right \vert = - \cos \beta l \cdot \sinh \beta l - \cosh \beta l \cdot \sin \beta l

Изобразим график f(z)=coszsinhz+coshzsinzf(z) = \cos z \cdot \sinh z + \cosh z \cdot \sin z:

import numpy as np import matplotlib.pyplot as plt maxX = 9.0 x = np.linspace(0.0, maxX, 1000) f = np.cos(x) * np.sinh(x) + np.cosh(x)*np.sin(x) plt.plot(x, 1e-3*f) plt.plot([0, maxX], [0, 0], 'k--') plt.ylim([-0.001, 0.001]) plt.xlim([0.0, maxX]) plt.scatter([2.355, 5.49, 8.633], [0.0,0.0,0.0], color='red') plt.title("Plot of $f(z) = \cos z \cdot \sinh z + \cosh z \cdot \sin z$")
<matplotlib.text.Text at 0x7fbe12a42950>
Image in a Jupyter notebook

Так как по логике вещей β>0\beta > 0 и l>0l > 0, то нас интересуют только положительные корни уравнения. Приближённо найдём его корни и после этого уточним численно:

firstThreeRoots = [2.3651, 5.49835, 8.6396] print (np.cos(firstThreeRoots) * np.sinh(firstThreeRoots) + np.cosh(firstThreeRoots)*np.sin(firstThreeRoots))
[ -6.09976979e-04 9.43321371e-02 -8.79784566e-01]
from scipy.optimize import brenth import math as m r1 = brenth(lambda Z: m.cos(Z) * m.sinh(Z) + m.cosh(Z)*m.sin(Z), 2.32, 2.3651) r2 = brenth(lambda Z: m.cos(Z) * m.sinh(Z) + m.cosh(Z)*m.sin(Z), 5.42, 5.49835) r3 = brenth(lambda Z: m.cos(Z) * m.sinh(Z) + m.cosh(Z)*m.sin(Z), 8.56, 8.6396) for el in r1, r2, r3: print("Function value at root {} = {}".format(el, np.cos(el) * np.sinh(el) + np.cosh(el)*np.sin(el)))
Function value at root 2.36502037243 = 1.33226762955e-15 Function value at root 5.497803919 = 4.26325641456e-14 Function value at root 8.6393798287 = -4.09272615798e-12

Выразим теперь нетривиальное решение этой системы уравнений через βl\beta l: C4=C,    C2=CcoshβlcosβlC_4 = C, \;\; C_2 = C \frac{\cosh \beta l}{\cos \beta l}

Таким образом, собственные функции задаются формулой: Y(x)=C(coshβlcosβlcosβx+coshβx)Y(x) = C \left ( \frac{\cosh \beta l}{\cos \beta l} \cdot \cos \beta x + \cosh \beta x \right)

Для простоты теперь будем считать, что стержень имеет единичную длину, то есть l=1l = 1. Найдём аналитическую формулу для квадрата L2L^2-нормы собственной функции этой задачи: 0l(coshβlcosβlcosβx+coshβx)2dx\int\limits_{0}^{l} \left (\frac{\cosh \beta l}{\cos \beta l} \cdot \cos \beta x + \cosh \beta x \right)^2 dx

import sympy as sp from sympy import integrate sp.init_printing() x, l = sp.symbols('x l') b = sp.Symbol('beta', positive=True) funcNorm = sp.integrate((sp.cosh(b*x)+ (sp.cosh(b*l)*sp.cos(b*x)/sp.cos(b*l)))**2, (x, 0, l))

Значит, квадрат L2L^2-нормы этой функции равен

funcNorm
lsin2(βl)cosh2(βl)2cos2(βl)l2sinh2(βl)+lcosh2(βl)+3sin(βl)cosh2(βl)2βcos(βl)+32βsinh(βl)cosh(βl)\frac{l \sin^{2}{\left (\beta l \right )} \cosh^{2}{\left (\beta l \right )}}{2 \cos^{2}{\left (\beta l \right )}} - \frac{l}{2} \sinh^{2}{\left (\beta l \right )} + l \cosh^{2}{\left (\beta l \right )} + \frac{3 \sin{\left (\beta l \right )} \cosh^{2}{\left (\beta l \right )}}{2 \beta \cos{\left (\beta l \right )}} + \frac{3}{2 \beta} \sinh{\left (\beta l \right )} \cosh{\left (\beta l \right )}

Построим нормированные собственные функции:

xs = np.linspace(0.0, 1.0, 1000) r = r1 C = 1.0/m.sqrt(funcNorm.subs([(b, r), (l, 1.0)])) f = C*(m.cosh(r) * np.cos(r*xs) /m.cos(r) + np.cosh(r*xs)) plt.plot(xs, f, label='${C1} \cdot \cos{{({b}x)}}+\ {C2} \cdot \cosh{{({b}x)}}$'.format(C1 = C*m.cosh(r)/m.cos(r), C2 = C, b = r)) plt.title('First Normalized Eigen Form') plt.xlabel('$x$') plt.ylabel('$y$') plt.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15))
<matplotlib.legend.Legend at 0x7fb7f9a48b10>
Image in a Jupyter notebook
r = r2 C = 1.0/m.sqrt(funcNorm.subs([(b, r), (l, 1.0)])) f = C*(m.cosh(r) * np.cos(r*xs) /m.cos(r) + np.cosh(r*xs)) plt.plot(xs, f, label='${C1} \cdot \cos{{({b}x)}} +\ {C2} \cdot \cosh{{({b}x)}}$'.format(C1 = C*m.cosh(r)/m.cos(r), C2 = C, b = r)) plt.title('Second Normalized Eigen Form') plt.xlabel('$x$') plt.ylabel('$y$') plt.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15))
<matplotlib.legend.Legend at 0x7fb7f8f4e090>
Image in a Jupyter notebook
r = r3 C = 1.0/m.sqrt(funcNorm.subs([(b, r), (l, 1.0)])) f = C*(m.cosh(r) * np.cos(r*xs) /m.cos(r) + np.cosh(r*xs)) plt.plot(xs, f, label='${C1} \cdot \cos{{({b}x)}} +\ {C2} \cdot \cosh{{({b}x)}}$'.format(C1 = C*m.cosh(r)/m.cos(r), C2 = C, b = r)) plt.title('Third Normalized Eigen Form') plt.xlabel('$x$') plt.ylabel('$y$') plt.legend(loc='upper center', bbox_to_anchor=(0.5, -0.15))
<matplotlib.legend.Legend at 0x7fb7f7d58850>
Image in a Jupyter notebook

Формы деформации в полиномиальном виде

Аппроксимируем первую, вторую и третью собственные формы с помощью следующих полиномов:

Y0(x)=k=04akxk,  Y1(x)=l=05alxl,  Y2(x)=m=06amxm,  Y3(x)=p=07apxpY_0 (x) = \sum_{k=0}^{4} a_k x^k, \; Y_1 (x) = \sum_{l=0}^{5} a_l x^l, \; Y_2 (x) = \sum_{m=0}^{6} a_m x^m, \; Y_3 (x) = \sum_{p=0}^{7} a_p x^p

Тогда получится система уравнений (напомним, что мы уже приняли l=1l=1):

(условия ортонормированности) Y0,Y0=Y1,Y1=Y2,Y2=Y3,Y3=1\langle Y_0, Y_0 \rangle = \langle Y_1, Y_1 \rangle = \langle Y_2, Y_2 \rangle = \langle Y_3, Y_3 \rangle = 1 Y0,Y1=Y0,Y2=Y0,Y3=Y1,Y2=Y1,Y3=Y2,Y3=0\langle Y_0, Y_1 \rangle = \langle Y_0, Y_2 \rangle = \langle Y_0, Y_3 \rangle =\langle Y_1, Y_2 \rangle = \langle Y_1, Y_3 \rangle = \langle Y_2, Y_3 \rangle = 0 (граничные условия для каждого из полиномов) Y0(0)=Y0(0)=Y0(1)=Y0(1)=0Y_0'(0) = Y_0'''(0) = Y_0''(1) = Y_0'''(1) = 0 Y1(0)=Y1(0)=Y1(1)=Y1(1)=0Y_1'(0) = Y_1'''(0) = Y_1''(1) = Y_1'''(1) = 0 Y2(0)=Y2(0)=Y2(1)=Y2(1)=0Y_2'(0) = Y_2'''(0) = Y_2''(1) = Y_2'''(1) = 0 Y3(0)=Y3(0)=Y3(1)=Y3(1)=0Y_3'(0) = Y_3'''(0) = Y_3''(1) = Y_3'''(1) = 0

Получается 26 уравнений и 26 неизвестных. Следующий код задаёт систему этих уравнений для последующей обработки в системе компьютерной алгебры.

a0, a1, a2, a3, a4 = sp.symbols('a_0 a_1 a_2 a_3 a_4') b0, b1, b2, b3, b4, b5 = sp.symbols('b_0 b_1 b_2 b_3 b_4 b_5') c0, c1, c2, c3, c4, c5, c6 = sp.symbols('c_0 c_1 c_2 c_3 c_4 c_5 c_6') d0, d1, d2, d3, d4, d5, d6, d7 = sp.symbols('d_0 d_1 d_2 d_3 d_4 d_5 d_6 d_7') Y0 = a0+a1*x+a2*x**2+a3*x**3+a4*x**4 Y1 = b0+b1*x+b2*x**2+b3*x**3+b4*x**4+b5*x**5 Y2 = c0+c1*x+c2*x**2+c3*x**3+c4*x**4+c5*x**5+c6*x**6 Y3 = d0+d1*x+d2*x**2+d3*x**3+d4*x**4+d5*x**5+d6*x**6+d7*x**7 equations = [sp.integrate(Y0*Y0, (x, 0, 1))-1, sp.integrate(Y1*Y1, (x, 0, 1))-1, sp.integrate(Y2*Y2, (x, 0, 1))-1, sp.integrate(Y3*Y3, (x, 0, 1))-1, sp.integrate(Y0*Y1, (x, 0, 1)), sp.integrate(Y0*Y2, (x, 0, 1)), sp.integrate(Y0*Y3, (x, 0, 1)), sp.integrate(Y1*Y2, (x, 0, 1)), sp.integrate(Y1*Y3, (x, 0, 1)), sp.integrate(Y2*Y3, (x, 0, 1)), sp.diff(Y0, x).subs([(x, 0)]), sp.diff(Y1, x).subs([(x, 0)]), sp.diff(Y2, x).subs([(x, 0)]), sp.diff(Y3, x).subs([(x, 0)]), sp.diff(Y0, x, x, x).subs([(x, 0)]), sp.diff(Y1, x, x, x).subs([(x, 0)]), sp.diff(Y2, x, x, x).subs([(x, 0)]), sp.diff(Y3, x, x, x).subs([(x, 0)]), sp.diff(Y0, x, x).subs([(x, 1)]), sp.diff(Y1, x, x).subs([(x, 1)]), sp.diff(Y2, x, x).subs([(x, 1)]), sp.diff(Y3, x, x).subs([(x, 1)]), sp.diff(Y0, x, x, x).subs([(x, 1)]), sp.diff(Y1, x, x, x).subs([(x, 1)]), sp.diff(Y2, x, x, x).subs([(x, 1)]), sp.diff(Y3, x, x, x).subs([(x, 1)]),]

Так выглядит система уравнений для решения этой задачи в виде F=0F=0:

sp.Matrix(equations)
[a02+a0a1+2a03a2+a0a32+2a05a4+a123+a1a22+2a15a3+a1a43+a225+a2a33+2a27a4+a327+a3a44+a4291b02+b0b1+2b03b2+b0b32+2b05b4+b0b53+b123+b1b22+2b15b3+b1b43+2b17b5+b225+b2b33+2b27b4+b2b54+b327+b3b44+2b39b5+b429+b4b55+b52111c02+c0c1+2c03c2+c0c32+2c05c4+c0c53+2c07c6+c123+c1c22+2c15c3+c1c43+2c17c5+c1c64+c225+c2c33+2c27c4+c2c54+2c29c6+c327+c3c44+2c39c5+c3c65+c429+c4c55+2c411c6+c5211+c5c66+c62131d02+d0d1+2d03d2+d0d32+2d05d4+d0d53+2d07d6+d0d74+d123+d1d22+2d15d3+d1d43+2d17d5+d1d64+2d19d7+d225+d2d33+2d27d4+d2d54+2d29d6+d2d75+d327+d3d44+2d39d5+d3d65+2d311d7+d429+d4d55+2d411d6+d4d76+d5211+d5d66+2d513d7+d6213+d6d77+d72151a0b0+a0b12+a0b23+a0b34+a0b45+a0b56+a1b02+a1b13+a1b24+a1b35+a1b46+a1b57+a2b03+a2b14+a2b25+a2b36+a2b47+a2b58+a3b04+a3b15+a3b26+a3b37+a3b48+a3b59+a4b05+a4b16+a4b27+a4b38+a4b49+a4b510a0c0+a0c12+a0c23+a0c34+a0c45+a0c56+a0c67+a1c02+a1c13+a1c24+a1c35+a1c46+a1c57+a1c68+a2c03+a2c14+a2c25+a2c36+a2c47+a2c58+a2c69+a3c04+a3c15+a3c26+a3c37+a3c48+a3c59+a3c610+a4c05+a4c16+a4c27+a4c38+a4c49+a4c510+a4c611a0d0+a0d12+a0d23+a0d34+a0d45+a0d56+a0d67+a0d78+a1d02+a1d13+a1d24+a1d35+a1d46+a1d57+a1d68+a1d79+a2d03+a2d14+a2d25+a2d36+a2d47+a2d58+a2d69+a2d710+a3d04+a3d15+a3d26+a3d37+a3d48+a3d59+a3d610+a3d711+a4d05+a4d16+a4d27+a4d38+a4d49+a4d510+a4d611+a4d712b0c0+b0c12+b0c23+b0c34+b0c45+b0c56+b0c67+b1c02+b1c13+b1c24+b1c35+b1c46+b1c57+b1c68+b2c03+b2c14+b2c25+b2c36+b2c47+b2c58+b2c69+b3c04+b3c15+b3c26+b3c37+b3c48+b3c59+b3c610+b4c05+b4c16+b4c27+b4c38+b4c49+b4c510+b4c611+b5c06+b5c17+b5c28+b5c39+b5c410+b5c511+b5c612b0d0+b0d12+b0d23+b0d34+b0d45+b0d56+b0d67+b0d78+b1d02+b1d13+b1d24+b1d35+b1d46+b1d57+b1d68+b1d79+b2d03+b2d14+b2d25+b2d36+b2d47+b2d58+b2d69+b2d710+b3d04+b3d15+b3d26+b3d37+b3d48+b3d59+b3d610+b3d711+b4d05+b4d16+b4d27+b4d38+b4d49+b4d510+b4d611+b4d712+b5d06+b5d17+b5d28+b5d39+b5d410+b5d511+b5d612+b5d713c0d0+c0d12+c0d23+c0d34+c0d45+c0d56+c0d67+c0d78+c1d02+c1d13+c1d24+c1d35+c1d46+c1d57+c1d68+c1d79+c2d03+c2d14+c2d25+c2d36+c2d47+c2d58+c2d69+c2d710+c3d04+c3d15+c3d26+c3d37+c3d48+c3d59+c3d610+c3d711+c4d05+c4d16+c4d27+c4d38+c4d49+c4d510+c4d611+c4d712+c5d06+c5d17+c5d28+c5d39+c5d410+c5d511+c5d612+c5d713+c6d07+c6d18+c6d29+c6d310+c6d411+c6d512+c6d613+c6d714a1b1c1d16a36b36c36d32a2+6a3+12a42b2+6b3+12b4+20b52c2+6c3+12c4+20c5+30c62d2+6d3+12d4+20d5+30d6+42d76a3+24a46b3+24b4+60b56c3+24c4+60c5+120c66d3+24d4+60d5+120d6+210d7]\left[\begin{matrix}a_{0}^{2} + a_{0} a_{1} + \frac{2 a_{0}}{3} a_{2} + \frac{a_{0} a_{3}}{2} + \frac{2 a_{0}}{5} a_{4} + \frac{a_{1}^{2}}{3} + \frac{a_{1} a_{2}}{2} + \frac{2 a_{1}}{5} a_{3} + \frac{a_{1} a_{4}}{3} + \frac{a_{2}^{2}}{5} + \frac{a_{2} a_{3}}{3} + \frac{2 a_{2}}{7} a_{4} + \frac{a_{3}^{2}}{7} + \frac{a_{3} a_{4}}{4} + \frac{a_{4}^{2}}{9} - 1\\b_{0}^{2} + b_{0} b_{1} + \frac{2 b_{0}}{3} b_{2} + \frac{b_{0} b_{3}}{2} + \frac{2 b_{0}}{5} b_{4} + \frac{b_{0} b_{5}}{3} + \frac{b_{1}^{2}}{3} + \frac{b_{1} b_{2}}{2} + \frac{2 b_{1}}{5} b_{3} + \frac{b_{1} b_{4}}{3} + \frac{2 b_{1}}{7} b_{5} + \frac{b_{2}^{2}}{5} + \frac{b_{2} b_{3}}{3} + \frac{2 b_{2}}{7} b_{4} + \frac{b_{2} b_{5}}{4} + \frac{b_{3}^{2}}{7} + \frac{b_{3} b_{4}}{4} + \frac{2 b_{3}}{9} b_{5} + \frac{b_{4}^{2}}{9} + \frac{b_{4} b_{5}}{5} + \frac{b_{5}^{2}}{11} - 1\\c_{0}^{2} + c_{0} c_{1} + \frac{2 c_{0}}{3} c_{2} + \frac{c_{0} c_{3}}{2} + \frac{2 c_{0}}{5} c_{4} + \frac{c_{0} c_{5}}{3} + \frac{2 c_{0}}{7} c_{6} + \frac{c_{1}^{2}}{3} + \frac{c_{1} c_{2}}{2} + \frac{2 c_{1}}{5} c_{3} + \frac{c_{1} c_{4}}{3} + \frac{2 c_{1}}{7} c_{5} + \frac{c_{1} c_{6}}{4} + \frac{c_{2}^{2}}{5} + \frac{c_{2} c_{3}}{3} + \frac{2 c_{2}}{7} c_{4} + \frac{c_{2} c_{5}}{4} + \frac{2 c_{2}}{9} c_{6} + \frac{c_{3}^{2}}{7} + \frac{c_{3} c_{4}}{4} + \frac{2 c_{3}}{9} c_{5} + \frac{c_{3} c_{6}}{5} + \frac{c_{4}^{2}}{9} + \frac{c_{4} c_{5}}{5} + \frac{2 c_{4}}{11} c_{6} + \frac{c_{5}^{2}}{11} + \frac{c_{5} c_{6}}{6} + \frac{c_{6}^{2}}{13} - 1\\d_{0}^{2} + d_{0} d_{1} + \frac{2 d_{0}}{3} d_{2} + \frac{d_{0} d_{3}}{2} + \frac{2 d_{0}}{5} d_{4} + \frac{d_{0} d_{5}}{3} + \frac{2 d_{0}}{7} d_{6} + \frac{d_{0} d_{7}}{4} + \frac{d_{1}^{2}}{3} + \frac{d_{1} d_{2}}{2} + \frac{2 d_{1}}{5} d_{3} + \frac{d_{1} d_{4}}{3} + \frac{2 d_{1}}{7} d_{5} + \frac{d_{1} d_{6}}{4} + \frac{2 d_{1}}{9} d_{7} + \frac{d_{2}^{2}}{5} + \frac{d_{2} d_{3}}{3} + \frac{2 d_{2}}{7} d_{4} + \frac{d_{2} d_{5}}{4} + \frac{2 d_{2}}{9} d_{6} + \frac{d_{2} d_{7}}{5} + \frac{d_{3}^{2}}{7} + \frac{d_{3} d_{4}}{4} + \frac{2 d_{3}}{9} d_{5} + \frac{d_{3} d_{6}}{5} + \frac{2 d_{3}}{11} d_{7} + \frac{d_{4}^{2}}{9} + \frac{d_{4} d_{5}}{5} + \frac{2 d_{4}}{11} d_{6} + \frac{d_{4} d_{7}}{6} + \frac{d_{5}^{2}}{11} + \frac{d_{5} d_{6}}{6} + \frac{2 d_{5}}{13} d_{7} + \frac{d_{6}^{2}}{13} + \frac{d_{6} d_{7}}{7} + \frac{d_{7}^{2}}{15} - 1\\a_{0} b_{0} + \frac{a_{0} b_{1}}{2} + \frac{a_{0} b_{2}}{3} + \frac{a_{0} b_{3}}{4} + \frac{a_{0} b_{4}}{5} + \frac{a_{0} b_{5}}{6} + \frac{a_{1} b_{0}}{2} + \frac{a_{1} b_{1}}{3} + \frac{a_{1} b_{2}}{4} + \frac{a_{1} b_{3}}{5} + \frac{a_{1} b_{4}}{6} + \frac{a_{1} b_{5}}{7} + \frac{a_{2} b_{0}}{3} + \frac{a_{2} b_{1}}{4} + \frac{a_{2} b_{2}}{5} + \frac{a_{2} b_{3}}{6} + \frac{a_{2} b_{4}}{7} + \frac{a_{2} b_{5}}{8} + \frac{a_{3} b_{0}}{4} + \frac{a_{3} b_{1}}{5} + \frac{a_{3} b_{2}}{6} + \frac{a_{3} b_{3}}{7} + \frac{a_{3} b_{4}}{8} + \frac{a_{3} b_{5}}{9} + \frac{a_{4} b_{0}}{5} + \frac{a_{4} b_{1}}{6} + \frac{a_{4} b_{2}}{7} + \frac{a_{4} b_{3}}{8} + \frac{a_{4} b_{4}}{9} + \frac{a_{4} b_{5}}{10}\\a_{0} c_{0} + \frac{a_{0} c_{1}}{2} + \frac{a_{0} c_{2}}{3} + \frac{a_{0} c_{3}}{4} + \frac{a_{0} c_{4}}{5} + \frac{a_{0} c_{5}}{6} + \frac{a_{0} c_{6}}{7} + \frac{a_{1} c_{0}}{2} + \frac{a_{1} c_{1}}{3} + \frac{a_{1} c_{2}}{4} + \frac{a_{1} c_{3}}{5} + \frac{a_{1} c_{4}}{6} + \frac{a_{1} c_{5}}{7} + \frac{a_{1} c_{6}}{8} + \frac{a_{2} c_{0}}{3} + \frac{a_{2} c_{1}}{4} + \frac{a_{2} c_{2}}{5} + \frac{a_{2} c_{3}}{6} + \frac{a_{2} c_{4}}{7} + \frac{a_{2} c_{5}}{8} + \frac{a_{2} c_{6}}{9} + \frac{a_{3} c_{0}}{4} + \frac{a_{3} c_{1}}{5} + \frac{a_{3} c_{2}}{6} + \frac{a_{3} c_{3}}{7} + \frac{a_{3} c_{4}}{8} + \frac{a_{3} c_{5}}{9} + \frac{a_{3} c_{6}}{10} + \frac{a_{4} c_{0}}{5} + \frac{a_{4} c_{1}}{6} + \frac{a_{4} c_{2}}{7} + \frac{a_{4} c_{3}}{8} + \frac{a_{4} c_{4}}{9} + \frac{a_{4} c_{5}}{10} + \frac{a_{4} c_{6}}{11}\\a_{0} d_{0} + \frac{a_{0} d_{1}}{2} + \frac{a_{0} d_{2}}{3} + \frac{a_{0} d_{3}}{4} + \frac{a_{0} d_{4}}{5} + \frac{a_{0} d_{5}}{6} + \frac{a_{0} d_{6}}{7} + \frac{a_{0} d_{7}}{8} + \frac{a_{1} d_{0}}{2} + \frac{a_{1} d_{1}}{3} + \frac{a_{1} d_{2}}{4} + \frac{a_{1} d_{3}}{5} + \frac{a_{1} d_{4}}{6} + \frac{a_{1} d_{5}}{7} + \frac{a_{1} d_{6}}{8} + \frac{a_{1} d_{7}}{9} + \frac{a_{2} d_{0}}{3} + \frac{a_{2} d_{1}}{4} + \frac{a_{2} d_{2}}{5} + \frac{a_{2} d_{3}}{6} + \frac{a_{2} d_{4}}{7} + \frac{a_{2} d_{5}}{8} + \frac{a_{2} d_{6}}{9} + \frac{a_{2} d_{7}}{10} + \frac{a_{3} d_{0}}{4} + \frac{a_{3} d_{1}}{5} + \frac{a_{3} d_{2}}{6} + \frac{a_{3} d_{3}}{7} + \frac{a_{3} d_{4}}{8} + \frac{a_{3} d_{5}}{9} + \frac{a_{3} d_{6}}{10} + \frac{a_{3} d_{7}}{11} + \frac{a_{4} d_{0}}{5} + \frac{a_{4} d_{1}}{6} + \frac{a_{4} d_{2}}{7} + \frac{a_{4} d_{3}}{8} + \frac{a_{4} d_{4}}{9} + \frac{a_{4} d_{5}}{10} + \frac{a_{4} d_{6}}{11} + \frac{a_{4} d_{7}}{12}\\b_{0} c_{0} + \frac{b_{0} c_{1}}{2} + \frac{b_{0} c_{2}}{3} + \frac{b_{0} c_{3}}{4} + \frac{b_{0} c_{4}}{5} + \frac{b_{0} c_{5}}{6} + \frac{b_{0} c_{6}}{7} + \frac{b_{1} c_{0}}{2} + \frac{b_{1} c_{1}}{3} + \frac{b_{1} c_{2}}{4} + \frac{b_{1} c_{3}}{5} + \frac{b_{1} c_{4}}{6} + \frac{b_{1} c_{5}}{7} + \frac{b_{1} c_{6}}{8} + \frac{b_{2} c_{0}}{3} + \frac{b_{2} c_{1}}{4} + \frac{b_{2} c_{2}}{5} + \frac{b_{2} c_{3}}{6} + \frac{b_{2} c_{4}}{7} + \frac{b_{2} c_{5}}{8} + \frac{b_{2} c_{6}}{9} + \frac{b_{3} c_{0}}{4} + \frac{b_{3} c_{1}}{5} + \frac{b_{3} c_{2}}{6} + \frac{b_{3} c_{3}}{7} + \frac{b_{3} c_{4}}{8} + \frac{b_{3} c_{5}}{9} + \frac{b_{3} c_{6}}{10} + \frac{b_{4} c_{0}}{5} + \frac{b_{4} c_{1}}{6} + \frac{b_{4} c_{2}}{7} + \frac{b_{4} c_{3}}{8} + \frac{b_{4} c_{4}}{9} + \frac{b_{4} c_{5}}{10} + \frac{b_{4} c_{6}}{11} + \frac{b_{5} c_{0}}{6} + \frac{b_{5} c_{1}}{7} + \frac{b_{5} c_{2}}{8} + \frac{b_{5} c_{3}}{9} + \frac{b_{5} c_{4}}{10} + \frac{b_{5} c_{5}}{11} + \frac{b_{5} c_{6}}{12}\\b_{0} d_{0} + \frac{b_{0} d_{1}}{2} + \frac{b_{0} d_{2}}{3} + \frac{b_{0} d_{3}}{4} + \frac{b_{0} d_{4}}{5} + \frac{b_{0} d_{5}}{6} + \frac{b_{0} d_{6}}{7} + \frac{b_{0} d_{7}}{8} + \frac{b_{1} d_{0}}{2} + \frac{b_{1} d_{1}}{3} + \frac{b_{1} d_{2}}{4} + \frac{b_{1} d_{3}}{5} + \frac{b_{1} d_{4}}{6} + \frac{b_{1} d_{5}}{7} + \frac{b_{1} d_{6}}{8} + \frac{b_{1} d_{7}}{9} + \frac{b_{2} d_{0}}{3} + \frac{b_{2} d_{1}}{4} + \frac{b_{2} d_{2}}{5} + \frac{b_{2} d_{3}}{6} + \frac{b_{2} d_{4}}{7} + \frac{b_{2} d_{5}}{8} + \frac{b_{2} d_{6}}{9} + \frac{b_{2} d_{7}}{10} + \frac{b_{3} d_{0}}{4} + \frac{b_{3} d_{1}}{5} + \frac{b_{3} d_{2}}{6} + \frac{b_{3} d_{3}}{7} + \frac{b_{3} d_{4}}{8} + \frac{b_{3} d_{5}}{9} + \frac{b_{3} d_{6}}{10} + \frac{b_{3} d_{7}}{11} + \frac{b_{4} d_{0}}{5} + \frac{b_{4} d_{1}}{6} + \frac{b_{4} d_{2}}{7} + \frac{b_{4} d_{3}}{8} + \frac{b_{4} d_{4}}{9} + \frac{b_{4} d_{5}}{10} + \frac{b_{4} d_{6}}{11} + \frac{b_{4} d_{7}}{12} + \frac{b_{5} d_{0}}{6} + \frac{b_{5} d_{1}}{7} + \frac{b_{5} d_{2}}{8} + \frac{b_{5} d_{3}}{9} + \frac{b_{5} d_{4}}{10} + \frac{b_{5} d_{5}}{11} + \frac{b_{5} d_{6}}{12} + \frac{b_{5} d_{7}}{13}\\c_{0} d_{0} + \frac{c_{0} d_{1}}{2} + \frac{c_{0} d_{2}}{3} + \frac{c_{0} d_{3}}{4} + \frac{c_{0} d_{4}}{5} + \frac{c_{0} d_{5}}{6} + \frac{c_{0} d_{6}}{7} + \frac{c_{0} d_{7}}{8} + \frac{c_{1} d_{0}}{2} + \frac{c_{1} d_{1}}{3} + \frac{c_{1} d_{2}}{4} + \frac{c_{1} d_{3}}{5} + \frac{c_{1} d_{4}}{6} + \frac{c_{1} d_{5}}{7} + \frac{c_{1} d_{6}}{8} + \frac{c_{1} d_{7}}{9} + \frac{c_{2} d_{0}}{3} + \frac{c_{2} d_{1}}{4} + \frac{c_{2} d_{2}}{5} + \frac{c_{2} d_{3}}{6} + \frac{c_{2} d_{4}}{7} + \frac{c_{2} d_{5}}{8} + \frac{c_{2} d_{6}}{9} + \frac{c_{2} d_{7}}{10} + \frac{c_{3} d_{0}}{4} + \frac{c_{3} d_{1}}{5} + \frac{c_{3} d_{2}}{6} + \frac{c_{3} d_{3}}{7} + \frac{c_{3} d_{4}}{8} + \frac{c_{3} d_{5}}{9} + \frac{c_{3} d_{6}}{10} + \frac{c_{3} d_{7}}{11} + \frac{c_{4} d_{0}}{5} + \frac{c_{4} d_{1}}{6} + \frac{c_{4} d_{2}}{7} + \frac{c_{4} d_{3}}{8} + \frac{c_{4} d_{4}}{9} + \frac{c_{4} d_{5}}{10} + \frac{c_{4} d_{6}}{11} + \frac{c_{4} d_{7}}{12} + \frac{c_{5} d_{0}}{6} + \frac{c_{5} d_{1}}{7} + \frac{c_{5} d_{2}}{8} + \frac{c_{5} d_{3}}{9} + \frac{c_{5} d_{4}}{10} + \frac{c_{5} d_{5}}{11} + \frac{c_{5} d_{6}}{12} + \frac{c_{5} d_{7}}{13} + \frac{c_{6} d_{0}}{7} + \frac{c_{6} d_{1}}{8} + \frac{c_{6} d_{2}}{9} + \frac{c_{6} d_{3}}{10} + \frac{c_{6} d_{4}}{11} + \frac{c_{6} d_{5}}{12} + \frac{c_{6} d_{6}}{13} + \frac{c_{6} d_{7}}{14}\\a_{1}\\b_{1}\\c_{1}\\d_{1}\\6 a_{3}\\6 b_{3}\\6 c_{3}\\6 d_{3}\\2 a_{2} + 6 a_{3} + 12 a_{4}\\2 b_{2} + 6 b_{3} + 12 b_{4} + 20 b_{5}\\2 c_{2} + 6 c_{3} + 12 c_{4} + 20 c_{5} + 30 c_{6}\\2 d_{2} + 6 d_{3} + 12 d_{4} + 20 d_{5} + 30 d_{6} + 42 d_{7}\\6 a_{3} + 24 a_{4}\\6 b_{3} + 24 b_{4} + 60 b_{5}\\6 c_{3} + 24 c_{4} + 60 c_{5} + 120 c_{6}\\6 d_{3} + 24 d_{4} + 60 d_{5} + 120 d_{6} + 210 d_{7}\end{matrix}\right]

Решим аналитически с использованием системы компьютерной алгебры:

solutions = sp.solve(equations, [a0, a1, a2, a3, a4, b0, b1, b2, b3, b4, b5, c0, c1, c2, c3, c4, c5, c6, d0, d1, d2, d3, d4, d5, d6, d7], dict='True')

Таким образом, есть

len(solutions)
1616

решений и они выглядят следующим образом:

solutions[0]
{a0:1,a1:0,a2:0,a3:0,a4:0,b0:263294411,b1:0,b2:563294274,b3:0,b4:563294548,b5:63294274,c0:29183207149748661121381031,c1:0,c2:31955320714974866747587354,c3:0,c4:2621153207149748661495174708,c5:144487320714974866747587354,c6:1683207149748662728421,d0:346560115931879979493115289670545081693,d1:0,d2:15415360815931879979493115289670545081693,d3:0,d4:185310069015931879979493115289670545081693,d5:394972082415931879979493115289670545081693,d6:309395814015931879979493115289670545081693,d7:31215931879979493115106167833}\left \{ a_{0} : -1, \quad a_{1} : 0, \quad a_{2} : 0, \quad a_{3} : 0, \quad a_{4} : 0, \quad b_{0} : - \frac{2 \sqrt{63294}}{411}, \quad b_{1} : 0, \quad b_{2} : \frac{5 \sqrt{63294}}{274}, \quad b_{3} : 0, \quad b_{4} : - \frac{5 \sqrt{63294}}{548}, \quad b_{5} : \frac{\sqrt{63294}}{274}, \quad c_{0} : - \frac{2918 \sqrt{320714974866}}{1121381031}, \quad c_{1} : 0, \quad c_{2} : \frac{31955 \sqrt{320714974866}}{747587354}, \quad c_{3} : 0, \quad c_{4} : - \frac{262115 \sqrt{320714974866}}{1495174708}, \quad c_{5} : \frac{144487 \sqrt{320714974866}}{747587354}, \quad c_{6} : - \frac{168 \sqrt{320714974866}}{2728421}, \quad d_{0} : - \frac{3465601 \sqrt{15931879979493115}}{289670545081693}, \quad d_{1} : 0, \quad d_{2} : \frac{154153608 \sqrt{15931879979493115}}{289670545081693}, \quad d_{3} : 0, \quad d_{4} : - \frac{1853100690 \sqrt{15931879979493115}}{289670545081693}, \quad d_{5} : \frac{3949720824 \sqrt{15931879979493115}}{289670545081693}, \quad d_{6} : - \frac{3093958140 \sqrt{15931879979493115}}{289670545081693}, \quad d_{7} : \frac{312 \sqrt{15931879979493115}}{106167833}\right \}

Количество решений (16 штук) может быть объяснено тем, что если (Y^1,Y^2,Y^3,Y^4)(\hat{Y}_1, \hat{Y}_2, \hat{Y}_3, \hat{Y}_4) — четвёрка решений, удовлетворяющая этой системе уравнений, то (±Y^1,±Y^2,±Y^3,±Y^4)(\pm\hat{Y}_1, \pm\hat{Y}_2, \pm\hat{Y}_3, \pm \hat{Y}_4) — тоже решение по очевидным соображениям. Поэтому можно взять любое из решений этой системы и получить любое другое, скорректировав знак.

sampleSolution = solutions[0]

Тогда получаем следующие аппроксимации собственных форм в полиномиальном виде:

  • первая

Y1res = Y1.subs(sampleSolution) Y1res
63294x5274563294548x4+563294274x2263294411\frac{\sqrt{63294} x^{5}}{274} - \frac{5 \sqrt{63294}}{548} x^{4} + \frac{5 \sqrt{63294}}{274} x^{2} - \frac{2 \sqrt{63294}}{411}
  • вторая

Y2res = Y2.subs(sampleSolution) Y2res
1683207149748662728421x6+144487320714974866747587354x52621153207149748661495174708x4+31955320714974866747587354x229183207149748661121381031- \frac{168 \sqrt{320714974866}}{2728421} x^{6} + \frac{144487 \sqrt{320714974866}}{747587354} x^{5} - \frac{262115 \sqrt{320714974866}}{1495174708} x^{4} + \frac{31955 \sqrt{320714974866}}{747587354} x^{2} - \frac{2918 \sqrt{320714974866}}{1121381031}
  • третья

Y3res = Y3.subs(sampleSolution) Y3res
31215931879979493115106167833x7309395814015931879979493115289670545081693x6+394972082415931879979493115289670545081693x5185310069015931879979493115289670545081693x4+15415360815931879979493115289670545081693x2346560115931879979493115289670545081693\frac{312 \sqrt{15931879979493115}}{106167833} x^{7} - \frac{3093958140 \sqrt{15931879979493115}}{289670545081693} x^{6} + \frac{3949720824 \sqrt{15931879979493115}}{289670545081693} x^{5} - \frac{1853100690 \sqrt{15931879979493115}}{289670545081693} x^{4} + \frac{154153608 \sqrt{15931879979493115}}{289670545081693} x^{2} - \frac{3465601 \sqrt{15931879979493115}}{289670545081693}

Теперь сравним, как выглядят численная аппроксимация и аналитически вычисленные собственные функции. Для первой собственной формы:

xs = np.linspace(0.0, 1.0, 1001) r = r1 C = 1.0/m.sqrt(funcNorm.subs([(b, r), (l, 1.0)])) f1Theoretical = C*(m.cosh(r) * np.cos(r*xs) /m.cos(r) + np.cosh(r*xs)) Y1Num = sp.lambdify(x, Y1res) f1Numerical = [Y1Num(xx) for xx in xs] plt.plot(xs, f1Theoretical, 'b', label='Analytical') plt.plot(xs, f1Numerical, 'r--', label='Numerical') plt.xlabel('$x$') plt.ylabel('$y$') plt.title('First Normalized Eigen Form: Analytical vs Numerical') plt.legend()
<matplotlib.legend.Legend at 0x7fb7f7bc58d0>
Image in a Jupyter notebook

Вторая собственная форма:

r = r2 C = 1.0/m.sqrt(funcNorm.subs([(b, r), (l, 1.0)])) f2Theoretical = C*(m.cosh(r) * np.cos(r*xs) /m.cos(r) + np.cosh(r*xs)) Y2Num = sp.lambdify(x, -1.0*Y2res) f2Numerical = [Y2Num(xx) for xx in xs] plt.plot(xs, f2Theoretical, 'b', label='Theoretical') plt.plot(xs, f2Numerical, 'r--', label='Numerical') plt.xlabel('$x$') plt.ylabel('$y$') plt.title('Second Normalized Eigen Form: Analytical vs Numerical') plt.legend()
<matplotlib.legend.Legend at 0x7fbe21c180d0>
Image in a Jupyter notebook

Третья собственная форма:

r = r3 C = 1.0/m.sqrt(funcNorm.subs([(b, r), (l, 1.0)])) f3Theoretical = C*(m.cosh(r) * np.cos(r*xs) /m.cos(r) + np.cosh(r*xs)) Y3Num = sp.lambdify(x, Y3res) f3Numerical = [Y3Num(xx) for xx in xs] plt.plot(xs, f3Theoretical, 'b', label='Analytical') plt.plot(xs, f3Numerical, 'r--', label='Numerical') plt.title('Third Normalized Eigen Form: Analytical vs Numerical') plt.xlabel('$x$') plt.ylabel('$y$') plt.legend()
<matplotlib.legend.Legend at 0x7fbe21b61c50>
Image in a Jupyter notebook

Вычислим среднеквадратичное отклонение, взяв в качестве узловых точек равномерную сетку на отрезке [0,1]\lbrack 0, 1 \rbrack с шагом 10310^{-3}:

  • для первой собственной формы он равен

np.sum(np.dot(f1Theoretical-f1Numerical,f1Theoretical-f1Numerical))/len(xs)
3.1677827161e053.1677827161e-05
  • для второй собственной формы он равен

np.sum(np.dot(f2Theoretical-f2Numerical,f2Theoretical-f2Numerical))/len(xs)
0.00105080205710.0010508020571
  • для третьей собственной формы он равен

np.sum(np.dot(f3Theoretical-f3Numerical,f3Theoretical-f3Numerical))/len(xs)
0.005274241592180.00527424159218

Данный метод является достаточно точным и подходит для быстрой аппроксимации собственных форм. Наибольшая ошибка квадрата среднеквадратичного отклонения достигается на третьей собственной форме и равна 0.005274241592180.00527424159218.

Устойчивость первых двух собственных форм

Рассмотрим следующее уравнение: EI4y(x,t)x4+p2y(x,t)x2+m2y(x,t)t2=0EI \frac{\partial^4 y(x, t)}{\partial x^4} + p \frac{\partial^2 y(x, t)}{\partial x^2} + m \frac{\partial^2 y(x, t)}{\partial t^2} = 0 вместе с граничными условиями y(x,t)xx=0=3y(x,t)x3x=0=2y(x,t)x2x=l=3y(x,t)x3x=l=0\left . \frac{\partial y(x, t)}{\partial x} \right \vert_{x = 0} = \left . \frac{\partial^3 y(x, t)}{\partial x^3} \right \vert_{x = 0} = \left . \frac{\partial^2 y(x, t)}{\partial x^2} \right \vert_{x = l} = \left . \frac{\partial^3 y(x, t)}{\partial x^3} \right \vert_{x = l} =0

Заменой x=lϕx = l \phi, t=mEIl2τt = \sqrt{\frac{m}{EI}} l^2 \tau мы переходим к безразмерной форме уравнения: 4y(ϕ,τ)ϕ4+b2y(ϕ,τ)ϕ2+2y(ϕ,τ)τ2=0,\frac{\partial^4 y(\phi, \tau)}{\partial \phi^4} + b \frac{\partial^2 y(\phi, \tau)}{\partial \phi^2} + \frac{\partial^2 y(\phi, \tau)}{\partial \tau^2} = 0, где b=pl2EIb = \frac{pl^2}{EI} и граничные условия принимают вид: y(ϕ,τ)ϕϕ=0=3y(ϕ,τ)ϕ3ϕ=0=2y(ϕ,τ)ϕ2ϕ=1=3y(ϕ,τ)ϕ3ϕ=1=0.\left . \frac{\partial y(\phi, \tau)}{\partial \phi} \right \vert_{\phi = 0} = \left . \frac{\partial^3 y(\phi, \tau)}{\partial \phi^3} \right \vert_{\phi = 0} = \left . \frac{\partial^2 y(\phi, \tau)}{\partial \phi^2} \right \vert_{\phi = 1} = \left . \frac{\partial^3 y(\phi, \tau)}{\partial \phi^3} \right \vert_{\phi = 1} =0.

Теперь в качестве решения подставим y(ϕ,τ)=X1(ϕ)T1(ϕ,τ)+X2(ϕ)T2(ϕ,τ)y(\phi, \tau) = X_1(\phi) T_1(\phi, \tau) + X_2(\phi) T_2 (\phi, \tau), где функции X1X_1, X2X_2, T1T_1 и T2T_2 связаны с первыми двумя собственными формами данного уравнения: k=12XkIVTk+bk=12XkIITk+k=12XkT¨k=0\sum_{k=1}^{2} X_k^{IV}T_k + b \sum_{k=1}^{2}X_k^{II} T_k + \sum_{k=1}^{2}X_k \ddot{T}_k = 0

Умножим скалярно это соотношение на X1X_1 и X2X_2: {T¨1+T1X1IV+bX1II,  X1+T2X2IV+bX2II,  X1=0T¨2+T1X1IV+bX1II,  X2+T2X2IV+bX2II,  X2=0\left \lbrace \begin{array}{ccc} \ddot{T}_1 + T_1 \cdot \langle X_1^{IV} + b X_1^{II} ,\; X_1 \rangle + T_2 \cdot \langle X_2^{IV} + b X_2^{II} ,\; X_1 \rangle & = &0\\ \ddot{T}_2 + T_1 \cdot \langle X_1^{IV} + b X_1^{II} ,\; X_2 \rangle + T_2 \cdot \langle X_2^{IV} + b X_2^{II} ,\; X_2 \rangle & = &0 \end{array} \right .

Примем aij=XjIV+bXjII,Xi  a_{ij} = \langle X_j^{IV}+b X_j^{II} , X_i \; \rangle. Тогда систему можно записать в виде: T¨=AT, \ddot{T} = A T, где T=(T1,T2)T = (T_1, T_2), а A=(a11a12a21a22)A = \left ( \begin{array}{cc} -a_{11} & -a_{12}\\ -a_{21} & -a_{22} \end{array} \right ). Вопросы асимптотического поведения (а, следовательно, и устойчивости) определяются решениями следующего характеристического уравнения: det(Aλ2I)=0\det{(A - \lambda^2 I)} = 0

Теперь уточним, чему же именно равны коэффициенты матрицы AA. Напомним, какому соотношению удовлетворяет собственная форма X(ϕ)X(\phi). В методе разделения переменных мы подставляем функцию вида y(ϕ,τ)=X(ϕ)T(τ)y(\phi, \tau) = X(\phi) T(\tau) в исходное уравнение и получаем XIVT+bXIIT+XT¨=0X^{IV} T + b X^{II} T + X \ddot{T} = 0 XIV(ϕ)+bXII(ϕ)X(ϕ)=T¨(τ)T(τ)=λ,\frac{X^{IV}(\phi) + b X^{II}(\phi)}{X(\phi)} = -\frac{\ddot{T}(\tau)}{T(\tau)} = \lambda,

откуда XIV(ϕ)+bXII(ϕ)=λX(ϕ)X^{IV}(\phi) + b X^{II}(\phi) = \lambda X(\phi). Но тогда это означает, что aij=XjIV+bXjII,Xi=λjXj,Xi=λjδija_{ij} = \langle X^{IV}_j + b X^{II}_j, X_i \rangle = \langle \lambda_j X_j, X_i \rangle = \lambda_j \delta_{ij}, где δij\delta_{ij} — дельта Кронекера, и на самом деле A=(λ100λ2).A = \left ( \begin{array}{cc} -\lambda_1 & 0\\ 0 & -\lambda_2 \end{array} \right ).

Система уравнений распадается на два независимых уравнения T¨1=λ1T1\ddot{T}_1 = -\lambda_1 T_1 и T¨2=λ2T2\ddot{T}_2 = -\lambda_2 T_2. Вопрос асимптотики решений для них очевиден: если λi>0\lambda_i > 0, то в решении получаются гармонические функции, которые ограничены при всех τR\tau \in \mathbb{R}; если же λi<0\lambda_i < 0, то решение есть суперпозиция cosh(λiτ)\cosh(\sqrt{-\lambda_i} \tau) и sinh(λiτ)\sinh(\sqrt{-\lambda_i} \tau), которая ограничена только в случае тривиального решения.

Всё, к чему теперь сводится задача — определить знак собственных чисел λ1\lambda_1 и λ2\lambda_2. Проанализируем какими могут быть корни характеристического уравнения в зависимости от λ\lambda: p4+bp2λ=0,p^4 + b p^2 - \lambda = 0, (p2+b2)2=λ+b24.\left (p^2 + \frac{b}{2} \right)^2 = \lambda + \frac{b^2}{4}. Из физических соображений b>0b > 0, поэтому можно выделить следующие случаи:

  • если λ<b24\lambda < -\frac{b^2}{4}, то корни имеют вид ±re±iθ\pm r e^{\pm i \theta}, r>0r>0, θ[0,π2]\theta \in \lbrack 0, \frac{\pi}{2} \rbrack

  • если b24<λ<0-\frac{b^2}{4} < \lambda < 0, то корни имеют вид ±iω\pm i \omega и ±iθ\pm i \theta, ω>0\omega > 0, θ>0\theta > 0

  • если λ>0\lambda > 0, то корни имеют вид ±iω\pm i\omega и ±α\pm \alpha, α>0\alpha > 0, β>0\beta > 0

Для устойчивости достаточно выяснить, возможно ли разрешить задачу на собственные числа при λ>0\lambda > 0. Если определитель системы для решения граничных условий всегда будет отличен от нуля, то это точно свидетельствует о том, что все собственные числа отрицательны, а значит система неустойчива.