CoCalc Shared Filestmp / dimakov / stability-domains.htmlOpen in CoCalc with one click!
Author: Evgenij Grines
Views : 11
Description: Jupyter html version of tmp/dimakov/stability-domains.ipynb
stability-domains
In [2]:
import sympy as sp

sp.init_printing()

o = sp.Symbol('omega', positive=True)
a = sp.Symbol('alpha', positive=True)
A = sp.Symbol('A', positive=True)
B = sp.Symbol('B', positive=True)
g = sp.Symbol('gamma', positive=True)
b = sp.Symbol('b', positive=True)
l = sp.Symbol('lambda')
z,x= sp.symbols('z x')
In [3]:
z**4 + b*z**2 - l
Out[3]:
bz2λ+z4b z^{2} - \lambda + z^{4}
In [4]:
C1, C2, C3, C4 = sp.symbols('C_1 C_2 C_3 C_4')

Формулы для выражения корней через коэффициенты:

Случай ±iω\pm i \omega, ±iα\pm i \alpha

In [5]:
sp.collect(sp.expand((z-sp.I*a)*(z+sp.I*a)*(z-sp.I*o)*(z+sp.I*o)), z)
Out[5]:
α2ω2+z4+z2(α2+ω2)\alpha^{2} \omega^{2} + z^{4} + z^{2} \left(\alpha^{2} + \omega^{2}\right)
In [6]:
gen_sol_1 = C1*sp.cos(a*x)+C2*sp.sin(a*x)+C3*sp.cos(o*x)+C4*sp.sin(o*x)
lhs1 = sp.Matrix([sp.diff(gen_sol_1, x).subs([(x, 0)]),
          sp.diff(gen_sol_1, x, x, x).subs([(x, 0)]),
          sp.diff(gen_sol_1, x, x).subs([(x, 1)]),
          sp.diff(gen_sol_1, x, x, x).subs([(x, 1)]),])
chpoly1 = sp.det((lhs1.jacobian([C1, C2, C3, C4])))
(chpoly1)
Out[6]:
α6ω3sin(α)cos(ω)+α5ω4sin(ω)cos(α)+α4ω5sin(α)cos(ω)α3ω6sin(ω)cos(α)- \alpha^{6} \omega^{3} \sin{\left (\alpha \right )} \cos{\left (\omega \right )} + \alpha^{5} \omega^{4} \sin{\left (\omega \right )} \cos{\left (\alpha \right )} + \alpha^{4} \omega^{5} \sin{\left (\alpha \right )} \cos{\left (\omega \right )} - \alpha^{3} \omega^{6} \sin{\left (\omega \right )} \cos{\left (\alpha \right )}
In [7]:
chpolyImaginaryPair = sp.lambdify([a, o], chpoly1)
In [8]:
print(sp.solve([a**2+o**2 - b,
          a**2*o**2+l],[a, o], dict=True)[7])
{omega: sqrt(b/2 - sqrt(b**2 + 4*lambda)/2), alpha: sqrt(b/2 + sqrt(b**2 + 4*lambda)/2)}

Эти формулы работают и дают положительные числа при b24<λ<0-\frac{b^2}{4} < \lambda < 0

Случай ±iα\pm i \alpha, ±γ\pm \gamma

In [9]:
sp.collect(sp.expand((z-g)*(z+g)*(z-sp.I*a)*(z+sp.I*a)), z)
Out[9]:
α2γ2+z4+z2(α2γ2)- \alpha^{2} \gamma^{2} + z^{4} + z^{2} \left(\alpha^{2} - \gamma^{2}\right)
In [10]:
gen_sol_2 = C1*sp.cos(a*x)+C2*sp.sin(a*x)+C3*sp.cosh(g*x)+C4*sp.sinh(g*x)
(gen_sol_2)
Out[10]:
C1cos(αx)+C2sin(αx)+C3cosh(γx)+C4sinh(γx)C_{1} \cos{\left (\alpha x \right )} + C_{2} \sin{\left (\alpha x \right )} + C_{3} \cosh{\left (\gamma x \right )} + C_{4} \sinh{\left (\gamma x \right )}
In [11]:
lhs2 = sp.Matrix([sp.diff(gen_sol_2, x, x).subs([(x, 1), (C2, 0), (C4, 0)]),
          sp.diff(gen_sol_2, x, x, x).subs([(x, 1)]),])


chpoly2 = sp.det((lhs2.jacobian([C1, C3])))
print(chpoly2)
-alpha**3*gamma**2*sin(alpha)*cosh(gamma) - alpha**2*gamma**3*cos(alpha)*sinh(gamma)
In [12]:
chpolyComplexReal = sp.lambdify((a, g), chpoly2)
In [13]:
sp.solve([a**2-g**2 - b,
          l - a**2*g**2], [a, g], dict=True)
Out[13]:
[{α:b212b2+4λ,γ:b212b2+4λ},{α:b212b2+4λ,γ:b212b2+4λ},{α:b212b2+4λ,γ:b212b2+4λ},{α:b212b2+4λ,γ:b212b2+4λ},{α:b2+12b2+4λ,γ:b2+12b2+4λ},{α:b2+12b2+4λ,γ:b2+12b2+4λ},{α:b2+12b2+4λ,γ:b2+12b2+4λ},{α:b2+12b2+4λ,γ:b2+12b2+4λ}]\left [ \left \{ \alpha : - \sqrt{\frac{b}{2} - \frac{1}{2} \sqrt{b^{2} + 4 \lambda}}, \quad \gamma : - \sqrt{- \frac{b}{2} - \frac{1}{2} \sqrt{b^{2} + 4 \lambda}}\right \}, \quad \left \{ \alpha : - \sqrt{\frac{b}{2} - \frac{1}{2} \sqrt{b^{2} + 4 \lambda}}, \quad \gamma : \sqrt{- \frac{b}{2} - \frac{1}{2} \sqrt{b^{2} + 4 \lambda}}\right \}, \quad \left \{ \alpha : \sqrt{\frac{b}{2} - \frac{1}{2} \sqrt{b^{2} + 4 \lambda}}, \quad \gamma : - \sqrt{- \frac{b}{2} - \frac{1}{2} \sqrt{b^{2} + 4 \lambda}}\right \}, \quad \left \{ \alpha : \sqrt{\frac{b}{2} - \frac{1}{2} \sqrt{b^{2} + 4 \lambda}}, \quad \gamma : \sqrt{- \frac{b}{2} - \frac{1}{2} \sqrt{b^{2} + 4 \lambda}}\right \}, \quad \left \{ \alpha : - \sqrt{\frac{b}{2} + \frac{1}{2} \sqrt{b^{2} + 4 \lambda}}, \quad \gamma : - \sqrt{- \frac{b}{2} + \frac{1}{2} \sqrt{b^{2} + 4 \lambda}}\right \}, \quad \left \{ \alpha : - \sqrt{\frac{b}{2} + \frac{1}{2} \sqrt{b^{2} + 4 \lambda}}, \quad \gamma : \sqrt{- \frac{b}{2} + \frac{1}{2} \sqrt{b^{2} + 4 \lambda}}\right \}, \quad \left \{ \alpha : \sqrt{\frac{b}{2} + \frac{1}{2} \sqrt{b^{2} + 4 \lambda}}, \quad \gamma : - \sqrt{- \frac{b}{2} + \frac{1}{2} \sqrt{b^{2} + 4 \lambda}}\right \}, \quad \left \{ \alpha : \sqrt{\frac{b}{2} + \frac{1}{2} \sqrt{b^{2} + 4 \lambda}}, \quad \gamma : \sqrt{- \frac{b}{2} + \frac{1}{2} \sqrt{b^{2} + 4 \lambda}}\right \}\right ]

Эти формулы работают при λ>0\lambda > 0

In [14]:
print(sp.solve([a**2-g**2 - b,
          l - a**2*g**2], [a, g], dict=True)[7])
{gamma: sqrt(-b/2 + sqrt(b**2 + 4*lambda)/2), alpha: sqrt(b/2 + sqrt(b**2 + 4*lambda)/2)}

Случай ±a±bi\pm a \pm b i

In [15]:
sp.collect(sp.expand((z-(A+sp.I*B))*(z-(A-sp.I*B))*(z-(-A+sp.I*B))*(z-(-A-sp.I*B))), z)
Out[15]:
A4+2A2B2+B4+z4+z2(2A2+2B2)A^{4} + 2 A^{2} B^{2} + B^{4} + z^{4} + z^{2} \left(- 2 A^{2} + 2 B^{2}\right)

При λ<b24\lambda <-\frac{b^2}{4}

In [16]:
gen_sol_3 = C1*sp.exp(A*x)*sp.cos(B*x)+C2*sp.exp(A*x)*sp.sin(B*x)+C3*sp.exp(-A*x)*sp.cos(B*x)+C4*sp.exp(-A*x)*sp.sin(B*x)
(gen_sol_3)
Out[16]:
C1eAxcos(Bx)+C2eAxsin(Bx)+C3eAxcos(Bx)+C4eAxsin(Bx)C_{1} e^{A x} \cos{\left (B x \right )} + C_{2} e^{A x} \sin{\left (B x \right )} + C_{3} e^{- A x} \cos{\left (B x \right )} + C_{4} e^{- A x} \sin{\left (B x \right )}
In [17]:
lhs3 = sp.Matrix([sp.diff(gen_sol_3, x).subs([(x, 0)]),
          sp.diff(gen_sol_3, x, x, x).subs([(x, 0)]),
          sp.diff(gen_sol_3, x, x).subs([(x, 1)]),
          sp.diff(gen_sol_3, x, x, x).subs([(x, 1)]),])
(lhs3)
Out[17]:
[AC1AC3+BC2+BC4A3C1A3C3+3A2BC2+3A2BC43AB2C1+3AB2C3B3C2B3C4A2C1eAcos(B)+A2C2eAsin(B)+A2C3eAcos(B)+A2C4eAsin(B)2ABC1eAsin(B)+2ABC2eAcos(B)+2C3eAABsin(B)2C4eAABcos(B)B2C1eAcos(B)B2C2eAsin(B)B2C3eAcos(B)B2C4eAsin(B)A3C1eAcos(B)+A3C2eAsin(B)A3C3eAcos(B)A3C4eAsin(B)3A2BC1eAsin(B)+3A2BC2eAcos(B)3C3eAA2Bsin(B)+3C4eAA2Bcos(B)3AB2C1eAcos(B)3AB2C2eAsin(B)+3C3eAAB2cos(B)+3C4eAAB2sin(B)+B3C1eAsin(B)B3C2eAcos(B)+B3C3eAsin(B)B3C4eAcos(B)]\left[\begin{matrix}A C_{1} - A C_{3} + B C_{2} + B C_{4}\\A^{3} C_{1} - A^{3} C_{3} + 3 A^{2} B C_{2} + 3 A^{2} B C_{4} - 3 A B^{2} C_{1} + 3 A B^{2} C_{3} - B^{3} C_{2} - B^{3} C_{4}\\A^{2} C_{1} e^{A} \cos{\left (B \right )} + A^{2} C_{2} e^{A} \sin{\left (B \right )} + \frac{A^{2} C_{3}}{e^{A}} \cos{\left (B \right )} + \frac{A^{2} C_{4}}{e^{A}} \sin{\left (B \right )} - 2 A B C_{1} e^{A} \sin{\left (B \right )} + 2 A B C_{2} e^{A} \cos{\left (B \right )} + \frac{2 C_{3}}{e^{A}} A B \sin{\left (B \right )} - \frac{2 C_{4}}{e^{A}} A B \cos{\left (B \right )} - B^{2} C_{1} e^{A} \cos{\left (B \right )} - B^{2} C_{2} e^{A} \sin{\left (B \right )} - \frac{B^{2} C_{3}}{e^{A}} \cos{\left (B \right )} - \frac{B^{2} C_{4}}{e^{A}} \sin{\left (B \right )}\\A^{3} C_{1} e^{A} \cos{\left (B \right )} + A^{3} C_{2} e^{A} \sin{\left (B \right )} - \frac{A^{3} C_{3}}{e^{A}} \cos{\left (B \right )} - \frac{A^{3} C_{4}}{e^{A}} \sin{\left (B \right )} - 3 A^{2} B C_{1} e^{A} \sin{\left (B \right )} + 3 A^{2} B C_{2} e^{A} \cos{\left (B \right )} - \frac{3 C_{3}}{e^{A}} A^{2} B \sin{\left (B \right )} + \frac{3 C_{4}}{e^{A}} A^{2} B \cos{\left (B \right )} - 3 A B^{2} C_{1} e^{A} \cos{\left (B \right )} - 3 A B^{2} C_{2} e^{A} \sin{\left (B \right )} + \frac{3 C_{3}}{e^{A}} A B^{2} \cos{\left (B \right )} + \frac{3 C_{4}}{e^{A}} A B^{2} \sin{\left (B \right )} + B^{3} C_{1} e^{A} \sin{\left (B \right )} - B^{3} C_{2} e^{A} \cos{\left (B \right )} + \frac{B^{3} C_{3}}{e^{A}} \sin{\left (B \right )} - \frac{B^{3} C_{4}}{e^{A}} \cos{\left (B \right )}\end{matrix}\right]
In [18]:
chpoly3 = sp.det((lhs3.jacobian([C1, C2, C3, C4])))
print(chpoly3)
-8*A**8*B*sin(B)*cos(B) - 2*A**7*B**2*exp(2*A)*sin(B)**2 - 2*A**7*B**2*exp(2*A)*cos(B)**2 + 2*A**7*B**2*exp(-2*A)*sin(B)**2 + 2*A**7*B**2*exp(-2*A)*cos(B)**2 - 24*A**6*B**3*sin(B)*cos(B) - 6*A**5*B**4*exp(2*A)*sin(B)**2 - 6*A**5*B**4*exp(2*A)*cos(B)**2 + 6*A**5*B**4*exp(-2*A)*sin(B)**2 + 6*A**5*B**4*exp(-2*A)*cos(B)**2 - 24*A**4*B**5*sin(B)*cos(B) - 6*A**3*B**6*exp(2*A)*sin(B)**2 - 6*A**3*B**6*exp(2*A)*cos(B)**2 + 6*A**3*B**6*exp(-2*A)*sin(B)**2 + 6*A**3*B**6*exp(-2*A)*cos(B)**2 - 8*A**2*B**7*sin(B)*cos(B) - 2*A*B**8*exp(2*A)*sin(B)**2 - 2*A*B**8*exp(2*A)*cos(B)**2 + 2*A*B**8*exp(-2*A)*sin(B)**2 + 2*A*B**8*exp(-2*A)*cos(B)**2
In [19]:
chpolyAllComplex = sp.lambdify((A, B), chpoly3)
In [20]:
sp.solve([2*B**2-2*A**2 - b,
          (A**2+B**2)**2+l], [A, B])
Out[20]:
[(12b2λ,b4λ2),(12b2λ,b4λ2),(12b2λ,b4λ2),(12b2λ,b4λ2),(12b+2λ,b4+λ2),(12b+2λ,b4+λ2),(12b+2λ,b4+λ2),(12b+2λ,b4+λ2)]\left [ \left ( - \frac{1}{2} \sqrt{- b - 2 \sqrt{- \lambda}}, \quad - \sqrt{\frac{b}{4} - \frac{\sqrt{- \lambda}}{2}}\right ), \quad \left ( - \frac{1}{2} \sqrt{- b - 2 \sqrt{- \lambda}}, \quad \sqrt{\frac{b}{4} - \frac{\sqrt{- \lambda}}{2}}\right ), \quad \left ( \frac{1}{2} \sqrt{- b - 2 \sqrt{- \lambda}}, \quad - \sqrt{\frac{b}{4} - \frac{\sqrt{- \lambda}}{2}}\right ), \quad \left ( \frac{1}{2} \sqrt{- b - 2 \sqrt{- \lambda}}, \quad \sqrt{\frac{b}{4} - \frac{\sqrt{- \lambda}}{2}}\right ), \quad \left ( - \frac{1}{2} \sqrt{- b + 2 \sqrt{- \lambda}}, \quad - \sqrt{\frac{b}{4} + \frac{\sqrt{- \lambda}}{2}}\right ), \quad \left ( - \frac{1}{2} \sqrt{- b + 2 \sqrt{- \lambda}}, \quad \sqrt{\frac{b}{4} + \frac{\sqrt{- \lambda}}{2}}\right ), \quad \left ( \frac{1}{2} \sqrt{- b + 2 \sqrt{- \lambda}}, \quad - \sqrt{\frac{b}{4} + \frac{\sqrt{- \lambda}}{2}}\right ), \quad \left ( \frac{1}{2} \sqrt{- b + 2 \sqrt{- \lambda}}, \quad \sqrt{\frac{b}{4} + \frac{\sqrt{- \lambda}}{2}}\right )\right ]

Формулы, описывающие удобное решение:

In [21]:
(sp.solve([2*B**2-2*A**2 - b,
          (A**2+B**2)**2+l], [A, B], dict=True)[7])
Out[21]:
{A:12b+2λ,B:b4+λ2}\left \{ A : \frac{1}{2} \sqrt{- b + 2 \sqrt{- \lambda}}, \quad B : \sqrt{\frac{b}{4} + \frac{\sqrt{- \lambda}}{2}}\right \}

Построение карты собственных значений

In [22]:
import numpy as np

N = 3001
M = 2001
bs = np.linspace(0.0 , 1.0, N)
ls = np.linspace(-0.3, 0.01, M)
In [23]:
caseImaginaryPairsDets = np.zeros((M, N))
caseComplexRealDets = np.zeros((M, N))
caseAllComplexDets = np.zeros((M, N))
In [24]:
for j, bb in enumerate(bs):
    for i, ll in enumerate(ls):
        caseImaginaryPairsDets[i, j] = np.nan
        caseComplexRealDets[i, j] = np.nan
        caseAllComplexDets[i, j] = np.nan
        if ll > 0:
            GG,AA = np.sqrt(-bb/2 + np.sqrt(bb**2 + 4*ll)/2), np.sqrt(bb/2 + np.sqrt(bb**2 + 4*ll)/2)
            caseComplexRealDets[i, j] = chpolyComplexReal(AA, GG)
        elif ll < - bb**2/4.0:
            BB, AA = np.sqrt(bb/4 + np.sqrt(-ll)/2), np.sqrt(-bb + 2*np.sqrt(-ll))/2
            caseAllComplexDets[i, j] = chpolyAllComplex(AA, BB)
        elif ((ll > -bb**2/4.0) and (ll < 0)):
            OO, AA = np.sqrt(bb/2 - np.sqrt(bb**2 + 4*ll)/2),  np.sqrt(bb/2 + np.sqrt(bb**2 + 4*ll)/2)
            caseImaginaryPairsDets[i, j] = chpolyImaginaryPair(AA,OO)
In [25]:
import matplotlib.pyplot as plt

plt.contour(bs, ls, caseAllComplexDets, [-0.001, 0.001], colors='red')
#plt.contourf(bs, ls, caseComplexRealDets, [-0.001, 0.001], colors='blue')
plt.contour(bs, ls, caseImaginaryPairsDets, [-0.001, 0.001], colors='green')
plt.xlabel('$b$')
plt.ylabel('$\lambda$')
#plt.plot([0.0, 1.0], [0.0, 0.0], 'k--')
Out[25]:
<matplotlib.text.Text at 0x7faf931af2d0>
In [ ]: