#!/usr/bin/env python
# coding: utf-8

# # Прикладные дифференциальные уравнения
# ## Занятие 10
# *Илья Щуров*
# 
# Факультет компьютерных наук, Прикладная математики и информатика, 2021-22 учебный год
# 
# [Страница курса](http://math-info.hse.ru/2021-22/Прикладные_дифференциальные_уравнения)

# In[1]:


import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt
import math
get_ipython().run_line_magic('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]))


# 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)


# In[8]:


g(e2)


# In[9]:


g(e1 + e2)


# In[10]:


g(e1) + g(e2)


# In[11]:


np.isclose(g(e1 + e2), g(e1) + g(e2)).all()


# In[12]:


g(3 * e1)


# In[13]:


np.isclose(3 * g(e1), g(3 * e1)).all()


# In[14]:


np.linalg.det(np.array([g(e1), g(e2)]))


# $$
# \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()


# ### Задача
# Найти $\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)


# In[19]:


np.linalg.det(np.array([g(e1), g(e2)]))


# In[20]:


g(e1)


# In[21]:


g(e2)


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

# 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, :]), '.')