{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Домашнее задание №3\n",
    "\n",
    "Это домашнее задание посвящено численному и математическому моделированию. \n",
    "\n",
    "Вашу программу следует разделить на следующие логические этапы.\n",
    "\n",
    "1. Чтение конфигурации из конфигурационного файла (физические параметры системы, коэффициенты уравнения, параметры методов решения задачи) и т.д.;\n",
    "2. Получение численного решения задачи при заданных параметрах;\n",
    "3. Визуализация полученного решения.\n",
    "\n",
    "Ваша задача заключается не только в программировании, но и в анализе самой задачи и её решений, которые могут иметь качественно разный характер при разных параметрах. Необходимо подготовить к моменту сдачи несколько конфигурационных файлов, демонстрирующих максимально разнообразие. Конфигурационные файлы разумно представить в виде JSON.\n",
    "\n",
    "Численное моделирование следует осуществлять на основе методов библиотеки `SciPy`, а не реализовывать численные методы самостоятельно. В целях визуализации допускается использовать любую библиотеку. На каждом графике должны быть подписаны оси (`ax.set_xlabel` и т.д.). Если в одних осях построены графики нескольких функций, то следует составить легенду (`ax.set_legend`).\n",
    "\n",
    "Желательно протестировать вашу программу. Если решение вашей задачи при каких-то параметрах очевидно, то следует проверить, что ваша программа способно его найти (в пределах погрешности). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## 1. Вычисление моментов инерции\n",
    "\n",
    "Пусть дано тело с плотностью $\\rho(x, y, z)$ в пределах области $D$. \n",
    "\n",
    "Найти:\n",
    "\n",
    "1. Его центр масс\n",
    "2. Тензор инерции\n",
    "3. Его главные оси\n",
    "\n",
    "И визуализировать это тело с его главными осями. \n",
    "\n",
    "Предусмотреть тела следующих форм: параллелепипед, сфера, цилиндр, конус и тор. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Задача двух тел\n",
    "\n",
    "Записать систему дифференциальных уравнений для задачи двух тел и реализовать её численное решение. Построить траектории этих тел (опционально можно визуализировать решение анимацией)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Тепловое равновесие \n",
    "\n",
    "Рассмотрим стержень длинны $l$, на концах которого поддерживаются температуры $T_0$ и $T_l$ соответственно. \n",
    "\n",
    "```{figure} /_static/lecture_specific/hw_3/heat.png\n",
    ":scale: 70%\n",
    "```\n",
    "\n",
    "Если направить ось $Ox$ вдоль длинны стержня как на рисунке выше и считать, что температура зависит только от координаты $x$ и что стержень теплоизолирован по его длине, то распределение температуры по длине стержня при достижении теплового равновесия описывается уравнением\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "\\dfrac{d}{dx}(\\rho(x)\\dfrac{dT}{dx}) = 0 \\\\\n",
    "T\\big|_{x=0} = T_0 \\\\\n",
    "T\\big|_{x=l} = T_l.\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Здесь $\\rho$ --- коэффициент теплопроводности. В общем случае этот коэффициент зависит от $x$, а в случае однородного стержня $\\rho(x)$ --- константа.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Осциллятор Ван дер Поля\n",
    "\n",
    "Реализовать численное решение уравнения Ван дер Поля\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "y'' - k (1- y^2)y' + y = 0 \\\\\n",
    "y\\big|_{t=0} = y_0 \\\\\n",
    "y'\\big|_{t=0} = y'_0\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Построить график решения в осях $(t, y)$ и $(y, y')$. Исследовать решение в зависимости от различных начальных условий, значений коэффициента $k$. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Уравнение Дуффинга\n",
    "\n",
    "В теории колебаний одним из наиболее известных нелинейных уравнений является уравнение Дуффинга, которое с учётом внешней силы записывается в виде\n",
    "\n",
    "$$\n",
    "\\ddot{x} + 2\\delta \\dot{x} + \\alpha x + \\beta x^3 = B \\sin(\\omega t)\n",
    "$$\n",
    "\n",
    "Пусть $\\delta=0.05$, $\\alpha=1$, $\\beta=1$, $\\omega=1$  \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Хищник-жертва\n",
    "\n",
    "Пусть $x$ --- плотность популяции жертвы, $y$ --- плотность популяции хищника, тогда, согласно модели Лотки-Вольтерра, модели хищник-жертва, можно записать следующие уравнения:\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "\\dot{x} = ax - bxy\\\\\n",
    "\\dot{y} = cxy - dy\n",
    "\\end{cases},\n",
    "$$\n",
    "\n",
    "где $a, b, c$ и $d$ положительные константы. Коэффициент $a$ описывает интенсивность размножения жертв, $b$ --- выедание хищниками жертв, $c$ --- увеличение биомассы хищников за счет выедания жертв, $d$ --- интенсивность естественной смерти хищников.\n",
    "\n",
    "Реализовать численное решение этой системы ОДУ и подобрать такие начальные плотности популяций и коэффициенты  $a, b, c$ и $d$, чтобы наблюдались периодические колебания, описывающие совместное проживание популяций жертв и хищников."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Аттрактор Лоренца\n",
    "\n",
    "Рассмотрим систему нелинейных ОДУ\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "\\dot{x} = \\sigma (y - x)\\\\\n",
    "\\dot{y} = x (r - z) - y\\\\\n",
    "\\dot{z} = xy - bz. \n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Реализовать численное решение этой системы, продемонстрировать странный аттрактор."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8. Заряды в электростатическом поле\n",
    "\n",
    "Пусть дано внешнее [электростатическое поле](https://ru.wikipedia.org/wiki/%D0%AD%D0%BB%D0%B5%D0%BA%D1%82%D1%80%D0%BE%D1%81%D1%82%D0%B0%D1%82%D0%B8%D1%87%D0%B5%D1%81%D0%BA%D0%B8%D0%B9_%D0%BF%D0%BE%D1%82%D0%B5%D0%BD%D1%86%D0%B8%D0%B0%D0%BB), потенциальная энергия взаимодействия которого с зарядом $q$ имеющим координаты $\\vec{r}=(x, y, z)$ имеет вид:\n",
    "\n",
    "$$\n",
    "W_q(\\vec{r}) = kq(x^2 + \\dfrac{y^2}{4} + \\dfrac{z^2}{16}).\n",
    "$$\n",
    "\n",
    "Предположим в это поле запускаются $N$ частиц с положительными зарядами $q_1, \\cdots, q_N$. Будем считать, что зависимость потенциальной энергии такой системы зарядов от положений этих зарядов $\\vec{r_1}, ..., \\vec{r_N}$ имеет вид\n",
    "\n",
    "$$\n",
    "W(\\overrightarrow{r_1}, ..., \\overrightarrow{r_N}) = \\sum_{i=1}^N W_{q_i}(\\overrightarrow{r_i}) + \\dfrac{1}{2}\\sum_{i\\neq j}k \\dfrac{q_i q_j}{r_{ij}},\n",
    "$$\n",
    "\n",
    "где $r_{ij}$ расстояние от $i$ до $j$ заряда. Т.е. потенциальная энергия системы является суммой энергии взаимодействия каждого заряда с электростатическим полем и энергией взаимодействия частиц между собой. \n",
    "\n",
    "Найти положение равновесия такой системы, считая, что оно достигается при минимальной потенциальной энергии: \n",
    "\n",
    "$$\n",
    "W(\\overrightarrow{r_1}, ..., \\overrightarrow{r_N}) \\sim \\min_{\\overrightarrow{r_1}, ..., \\overrightarrow{r_N}}.\n",
    "$$\n",
    "\n",
    "Визуализировать частицы в найденном положении равновесия."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 9. Двойной маятник\n",
    "\n",
    "Смоделируйте поведение двойного маятника.\n",
    "\n",
    "```{figure} /_static/lecture_specific/hw_3/dbl_pendulum.png\n",
    ":scale: 70%\n",
    "```\n",
    "\n",
    "Состояние системы в любой момент времени может быть описано углами $\\theta_1$ и $\\theta_2$ и угловыми скоростями $\\omega_1 = \\dot{\\theta_1}$ и $\\omega_2 = \\dot{\\theta_2}$. Если выбрать  углы в соответствии с рисунком выше, то уравнения движения имеют вид\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "\\varepsilon_1 = \\dfrac{-g(2m_1 + m_2)\\sin\\theta_1-gm_2\\sin(\\theta_1 - 2\\theta_2) - 2\\sin(\\theta_1 - \\theta_2)m_2(\\omega_2^2L_2 + \\omega_1^2l_1\\cos(\\theta_1 - \\theta_2))}{L_1(2m_1 + m_2 - m_2 \\cos(2\\theta_1 - 2 \\theta_2))}, \n",
    "\\\\\n",
    "\\varepsilon_2 = \\dfrac{2\\sin(\\theta_1 - \\theta_2)(\\omega_1^2L_1(m_1 + m_2) + g(m_1 + m_2)\\cos\\theta_1 + \\omega_2^2 L_2 m_2 \\cos(\\theta_1 - \\theta_2))}{L_1(2m_1 + m_2 - m_2 \\cos(2\\theta_1 - 2 \\theta_2))}.\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Здесь $\\varepsilon_1 = \\ddot{\\theta}_1$ и $\\varepsilon_2 = \\ddot{\\theta}_2$ --- угловые ускорения."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 10. Два маятника на пружине\n",
    "\n",
    "Рассмотрим механическую систему, состоящую из двух пружинных маятников, соединенных между собой ещё одной пружиной. Предположим, что в состоянии покоя системы все три пружины не испытывают деформаций.\n",
    "\n",
    "```{figure} /_static/lecture_specific/hw_3/springs.png\n",
    ":scale: 90%\n",
    "```\n",
    "\n",
    "Теперь предположим, что эту систему вывели из состояния покоя, сместив грузики из положений равновесий и придав им некоторые скорости. Пусть $x_1$ --- смещение первого грузика вправо относительно его положения равновесия, $x_2$ --- смещение второго грузика вправо относительно его положения равновесия, тогда уравнения движения имеют вид\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "ma_1 = -k_1 x_1 + k(x_2 - x_1),\\\\\n",
    "ma_2 = -k_2 x_2 - k(x_2 - x_1),\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "где $a_1 = \\ddot{x}_1$ и $a_2 = \\ddot{x}_2$ --- ускорения грузиков. Смоделировать поведение такой системы. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 11. Заряд в электромагнитном поле\n",
    "\n",
    "\n",
    "Точечный заряд движется в постоянном электромагнитном поле, т.е. него действует [сила Лоренца](https://ru.wikipedia.org/wiki/%D0%A1%D0%B8%D0%BB%D0%B0_%D0%9B%D0%BE%D1%80%D0%B5%D0%BD%D1%86%D0%B0#:~:text=%D0%A1%D0%B8%D0%BB%D0%B0%20%D0%9B%D0%BE%D1%80%D0%B5%D0%BD%D1%86%D0%B0%20%E2%80%94%20%D1%81%D0%B8%D0%BB%D0%B0%2C%20%D1%81,%D1%81%D0%B8%D1%81%D1%82%D0%B5%D0%BC%D0%B5%20%D0%B5%D0%B4%D0%B8%D0%BD%D0%B8%D1%86%20(%D0%A1%D0%98)%20%D0%B2%D1%8B%D1%80%D0%B0%D0%B6%D0%B0%D0%B5%D1%82%D1%81%D1%8F%20%D0%BA%D0%B0%D0%BA%3A)\n",
    "\n",
    "$$\n",
    "\\overrightarrow{F} = q(\\overrightarrow{E} + [\\overrightarrow{v} \\times \\overrightarrow{B}]).\n",
    "$$\n",
    "\n",
    "Здесь $q$ --- величина заряда, $\\overrightarrow{E}=\\overrightarrow{E}(\\vec{r})$ и $\\overrightarrow{B}=\\overrightarrow{B}(r)$ --- напряженности электрического и магнитного полей, $\\vec{r}$ и $\\vec{v}$ --- радиус-вектор положения частицы и вектор мгновенной скорости частицы.\n",
    "\n",
    "Пренебрегая силой тяжести, смоделировать полет заряда в таком поле решив уравнения движения\n",
    "\n",
    "$$\n",
    "m\\vec{a} = \\overrightarrow{F},\n",
    "$$\n",
    "при некоторых начальных условиях $\\vec{r}\\big|_{t=0} = \\vec{r}_0, \\, \\vec{v}\\big|_{t=0} = \\vec{v}_0$. Предусмотреть равномерные ($\\overrightarrow{E}(\\vec{r}) = \\overrightarrow{E}_0$, $\\overrightarrow{B}(\\vec{r}) = \\overrightarrow{B}_0$) и радиальные ($\\overrightarrow{E}(\\vec{r}) = \\vec{r}||\\overrightarrow{E}||/||\\vec{r}||$, $\\overrightarrow{B}(\\vec{r}) = \\vec{r}||\\overrightarrow{B}||/||\\vec{r}||$) электрические и магнитные поля."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12. Ядерный реактор\n",
    "\n",
    "Динамика ядерного реактора в простейшей модели описывается системой уравнений на мощность $N$ и температуру $T$ реактора:\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "\\dfrac{dN}{dt} = - \\dfrac{\\alpha (T - T_1)}{\\tau}N, \\\\\n",
    "C \\dfrac{dT}{dt} = N - k (T-T_0),\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "где $\\tau$ --- время жизни одного поколения нейтронов, $\\alpha > 0$ --- реактивность реактора, $C$ --- теплоемкость реактора, $T_1$ --- температура окружающей среды, $T_0 < T_1$ --- температура охладителя реактора, $k$ --- коэффициент теплоотдачи реактора.\n",
    "\n",
    "## 13. Химический реактор\n",
    "Реактор представляет собой открытую систему полного перемешивания с непрерывным тепло- и массообменом с окружающей средой. Химическая реакция протекает гомогенно с выделением тепла, которое отводится как потоком вещества, так и теплопередачей через стенки реактора. Зависимость скорости реакции от концентрации реагента и температуры определяется законом действующих масс Аррениуса. В этой простейшей модели динамика реактора описывается системой уравнений на концентрацию реагента $x$ и его температуру $y$:\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "\\dfrac{dx}{dt} = -x e^{-1/y} + \\lambda (x_0 - x), \\\\\n",
    "\\dfrac{dy}{dt} = -x e^{-1/y} + \\mu (y_0 - y),\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "где $x_0$ –-- концентрация подающегося реагента, $y_0$ – равновесная температура реактора при потушенной реакции, $\\lambda$ – скорость вывода реагента из реактора, $\\mu$ – скорость охлаждения за счёт вывода реагента и отвода тепла. \n"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}