1. cs-ode-appl-2021
  2. Lesson10.ipynb

Прикладные дифференциальные уравнения

Занятие 10

Илья Щуров

Факультет компьютерных наук, Прикладная математики и информатика, 2021-22 учебный год

Страница курса

In [1]:
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math
%matplotlib inline

Осциллятор с переменной жёсткостью

$$\ddot x = - (\sin t + 2) x$$$$\dot x = y, \quad \dot y = -(\sin t + 2) x$$
In [2]:
def osc(X, t):
    x, y = X
    return np.array([y, (-math.sin(2 * t) + 2) * x])
In [3]:
odeint(osc, np.array([1, 2]), np.array([0, 2]))
Out[3]:
array([[ 1.        ,  2.        ],
       [14.64673206, 21.63805004]])
In [4]:
def phase_flow(system, t_0, t_1):
    def g(x):
        return odeint(system, x, np.array([t_0, t_1]))[-1]
    return g
In [5]:
g = phase_flow(osc, 1, 4)

Линейность отображения фазового потока

In [6]:
e1 = np.array([1, 0])
e2 = np.array([0, 1])
In [7]:
g(e1)
Out[7]:
array([36.77682016, 42.06341281])
In [8]:
g(e2)
Out[8]:
array([29.09520864, 33.3047762 ])
In [9]:
g(e1 + e2)
Out[9]:
array([65.87203428, 75.36819653])
In [10]:
g(e1) + g(e2)
Out[10]:
array([65.8720288 , 75.36818901])
In [11]:
np.isclose(g(e1 + e2), g(e1) + g(e2)).all()
Out[11]:
True
In [12]:
g(3 * e1)
Out[12]:
array([110.33046451, 126.19024232])
In [13]:
np.isclose(3 * g(e1), g(3 * e1)).all()
Out[13]:
True
In [14]:
np.linalg.det(np.array([g(e1), g(e2)]))
Out[14]:
0.9999930497444511
$$ \dot x = tx + 2y, \quad \dot y = x \sin t - 3y $$
In [15]:
def linear_system(X, t):
    x, y = X
    return np.array([x * t + 2 * y,
              x * math.sin(t) - 3 * y])
In [16]:
g = phase_flow(linear_system, 0, 3)
In [17]:
np.isclose(g(e1 + e2), g(e1) + g(e2)).all()
Out[17]:
True

Задача

Найти $\det g_0^3$ аналитически.

$$\mathop{\mathrm{Tr}} A = t - 3$$$$\int_0^3 (t-3) dt = 9/2 - 9=-9/2=-4.5$$$$\det g_0^3 = e^{-4.5}$$
In [18]:
math.exp(3 ** 2 / 2 - 3 * 3)
Out[18]:
0.011108996538242306
In [19]:
np.linalg.det(np.array([g(e1), g(e2)]))
Out[19]:
0.011194714438770365
In [20]:
g(e1)
Out[20]:
array([196.84929067,   9.89932211])
In [21]:
g(e2)
Out[21]:
array([115.89484579,   5.82827401])

Сохранение фазового объёма

In [22]:
def pendulum(X, t):
    x, y = X
    return np.array([y,
                     -math.sin(x)])
In [23]:
initial_points = [x for x in np.random.uniform(-3, 3, size=(100000, 2))
                  if np.linalg.norm(x - np.array([0.8, 0.7])) < 0.4]
In [24]:
T = np.linspace(0, 50)
trajectories = [odeint(pendulum, initial_point, T) 
                for initial_point in initial_points]
In [25]:
from ipywidgets import interact
In [26]:
@interact(t=(0, 50), manual=True)
def draw_phase_volume(t):
    X = np.linspace(-8, 8, 500)
    Y = np.linspace(-4, 4, 500)
    x, y = np.meshgrid(X, Y)
    plt.contour(x, y, y ** 2 / 2 - np.cos(x), levels=20)
    plt.plot(*zip(*np.array(trajectories)[:, t, :]), '.')