CoCalc Public Filestmp / dimakov / dimakov-stability.ipynbOpen with one click!
Author: Evgenij Grines
Views : 54
Description: Jupyter notebook tmp/dimakov/dimakov-stability.ipynb
Compute Environment: Ubuntu 18.04 (Deprecated)

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

Рассмотрим следующее уравнение: 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:

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

Численное определение области устойчивости

В качестве аппроксимации области устойчивости воспользуемся следующим методом. Вместо настоящих первых двух мод будем брать их полиномиальные аппроксимации X^i\hat{X}_i, а матрица AA будет состоять из элементов aij=X^jIV+bX^jII,X^ia_{ij} = \langle \hat{X}^{IV}_j + b \hat{X}^{II}_j, \hat{X}_i \rangle. Зная собственные числа матрицы AA для системы T¨=AT\ddot{T}=AT, мы можем ответить на вопрос об устойчивости. А именно, если матрица AA имеет положительное собственное число pp, то есть и положительное собственное число λ=p\lambda = \sqrt{p} в системе дифференциальных уравнений T¨=AT\ddot{T} = AT, что приведёт к неустойчивости решений.

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

In [1]:
import sympy as sp from sympy import integrate sp.init_printing() x = sp.symbols('x') b = sp.symbols('b') 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') 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 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(Y0*Y1, (x, 0, 1)), sp.integrate(Y0*Y2, (x, 0, 1)), sp.integrate(Y1*Y2, (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(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(Y0, x, x).subs([(x, 1)]), sp.diff(Y1, x, x).subs([(x, 1)]), sp.diff(Y2, 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)]),] solutions = sp.solve(equations, [a0, a1, a2, a3, a4, b0, b1, b2, b3, b4, b5, c0, c1, c2, c3, c4, c5, c6,], dict='True') sampleSolution = solutions[0] Y1res = Y1.subs(sampleSolution) Y2res = Y2.subs(sampleSolution)

Так как первое решение вообще почти никакого отношения к реальным модам не имеет, то мы его можем отбросить и вычислить производные только для X^1\hat{X}_1 и X^2\hat{X}_2:

In [2]:
a11 = sp.integrate((sp.diff(Y1res, x, 4)+b*sp.diff(Y1res, x, 2))*Y1res ,(x, 0, 1)) a12 = sp.integrate((sp.diff(Y2res, x, 4)+b*sp.diff(Y2res, x, 2))*Y1res ,(x, 0, 1)) a21 = sp.integrate((sp.diff(Y1res, x, 4)+b*sp.diff(Y1res, x, 2))*Y2res ,(x, 0, 1)) a22 = sp.integrate((sp.diff(Y2res, x, 4)+b*sp.diff(Y2res, x, 2))*Y2res ,(x, 0, 1))

В итоге это приведёт к следующей аппроксимации матрицы AA (которая сейчас зависит от параметра bb):

In [3]:
A = sp.Matrix([[-a11, -a12],[-a21, -a22]]) (A.evalf())

Изобразим на графике след и определитель матрицы AA в зависимости от параметра bb:

In [7]:
import matplotlib.pyplot as plt import numpy as np bs = np.linspace(0, 12, 5000) dets = [sp.det(A).subs([(b, bb)]).evalf() for bb in bs] trfunc = A[0, 0]+A[1, 1] trs = [trfunc.subs([(b, bb)]).evalf() for bb in bs] plt.plot(bs, dets, label='det') plt.plot(bs, trs, label='tr') plt.xlabel('b') plt.plot([0.0, 12.0],[0.0, 0.0],'k--') plt.legend()
<matplotlib.legend.Legend at 0x7f26487ce210>

Мы наблюдаем, что до некоторого критического значения bb собственные числа матрицы A(b)A(b) либо одного знака, либо комплексные, но их сумма имеет тот же знак, что и след матрицы A(b)A(b), а он отрицателен. Это гарантирует отрицательность вещественных частей собственных чисел матрицы A(b)A(b).

In [8]:
import matplotlib.pyplot as plt import numpy as np bs = np.linspace(9.9, 9.95, 1000) dets = [sp.det(A).subs([(b, bb)]).evalf() for bb in bs] trfunc = A[0, 0]+A[1, 1] trs = [trfunc.subs([(b, bb)]).evalf() for bb in bs] plt.plot(bs, dets, label='det') #plt.plot(bs, trs, label='tr') #discr = np.array(trs)*np.array(trs) - 4.0*np.array(dets) #plt.plot(bs, discr, label='discr') plt.xlabel('b') plt.plot([9.9, 9.95],[0.0, 0.0],'k--') plt.legend()
<matplotlib.legend.Legend at 0x7f26440cf1d0>

Смена устойчивости происходит при b9.92b \approx 9.92, после чего определитель меняет знак и собственные числа точно разных знаков, что приводит к неустойчивости.

Однако, есть ещё одна смена устойчивости при b40.6b \approx 40.6 — после этого значения двухмодовое решение снова становится устойчивым:

In [9]:
import matplotlib.pyplot as plt import numpy as np bs = np.linspace(40.6, 40.8, 5000) dets = [sp.det(A).subs([(b, bb)]).evalf() for bb in bs] trfunc = A[0, 0]+A[1, 1] trs = [trfunc.subs([(b, bb)]).evalf() for bb in bs] plt.plot(bs, dets, label='det') plt.plot(bs, trs, label='tr') plt.xlabel('b') plt.plot([40.6, 40.8],[0.0, 0.0],'k--') plt.xlim([40.6, 40.8]) plt.legend()
<matplotlib.legend.Legend at 0x7f2643cc9590>
In [ ]:
In [ ]: