Open in CoCalc with one click!

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

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

  • 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) = 0

EImYIV(x)Y(x)=T(t)T=λ\frac{EI}{m} \cdot \frac{Y^{IV}(x)}{Y(x)} = -\frac{T''(t)}{T} = \lambda

EIYIV(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:

In [1]:
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>

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

In [2]:
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]
In [3]:
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

In [4]:
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-нормы этой функции равен

In [5]:
funcNorm

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

In [6]:
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>
In [7]:
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>
In [8]:
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>

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

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

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 неизвестных. Следующий код задаёт систему этих уравнений для последующей обработки в системе компьютерной алгебры.

In [9]:
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:

In [10]:
sp.Matrix(equations)

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

In [11]:
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')

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

In [12]:
len(solutions)

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

In [13]:
solutions[0]

Количество решений (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) — тоже решение по очевидным соображениям. Поэтому можно взять любое из решений этой системы и получить любое другое, скорректировав знак.

In [14]:
sampleSolution = solutions[0]

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

  • первая
In [15]:
Y1res = Y1.subs(sampleSolution) Y1res
  • вторая
In [16]:
Y2res = Y2.subs(sampleSolution) Y2res
  • третья
In [17]:
Y3res = Y3.subs(sampleSolution) Y3res

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

In [18]:
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>

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

In [19]:
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>

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

In [20]:
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>

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

  • для первой собственной формы он равен
In [21]:
np.sum(np.dot(f1Theoretical-f1Numerical,f1Theoretical-f1Numerical))/len(xs)
  • для второй собственной формы он равен
In [22]:
np.sum(np.dot(f2Theoretical-f2Numerical,f2Theoretical-f2Numerical))/len(xs)
  • для третьей собственной формы он равен
In [23]:
np.sum(np.dot(f3Theoretical-f3Numerical,f3Theoretical-f3Numerical))/len(xs)

Данный метод является достаточно точным и подходит для быстрой аппроксимации собственных форм. Наибольшая ошибка квадрата среднеквадратичного отклонения достигается на третьей собственной форме и равна 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. Если определитель системы для решения граничных условий всегда будет отличен от нуля, то это точно свидетельствует о том, что все собственные числа отрицательны, а значит система неустойчива.

In [ ]: