IntervalExpressions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Зависимость интервальных оценок областей значений выражений от их вида\n",
    "\n",
    "Результат вычисления арифметического выражения зависит от вида выражения, а именно от количества вхождений в него переменных. Эта особенность интервальных вычислений подробно обсуждается в первой главе книги [[1]](#SharyBook).\n",
    "\n",
    "Для иллюстрации приведём два практических примера интервальных расчётов для формул, выражающих естественнонаучные законы."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Пример 1. Уравнение Михаэлиса-Ментен \n",
    "\n",
    "Ферменты (называемые также энзимами) играют огромную роль в процессах обмена веществ в живых организмах, намного порядков повышая скорость химических реакций. Предложенное в 1913 году [уравнение Михаэлиса-Ментен](https://ru.wikipedia.org/wiki/%D0%A3%D1%80%D0%B0%D0%B2%D0%BD%D0%B5%D0%BD%D0%B8%D0%B5_%D0%9C%D0%B8%D1%85%D0%B0%D1%8D%D0%BB%D0%B8%D1%81%D0%B0_%E2%80%94_%D0%9C%D0%B5%D0%BD%D1%82%D0%B5%D0%BD) описывает зависимость скорости реакции, катализируемой ферментом, от концентрации субстрата:\n",
    "\n",
    "<a id='MichaelisMenten'></a>\n",
    "$$\\mathrm{v} = V_{max} \\frac{S}{S+K_M}, \\hspace{3 cm}(1)$$\n",
    "\n",
    "где $\\mathrm{v}$ --- скорость реакции, $V_{max} $ --- максимальная скорость реакции, \n",
    "$K_M$ --- константа Михаэлиса, $S$ ---  концентрация субстрата. \n",
    "\n",
    "Cкорость реакции $\\mathrm{v}$ будем выражать в единицах  $V_{max}$, поэтому зададимся значением $V_{max} =  1.$  \n",
    "\n",
    "Для определённости возьмём реакцию гидролиза аденозинтрифосфата, катализируемого \n",
    "миозином $K_M = 1.44 \\cdot 10^{-4} $ моль/л. \n",
    "\n",
    "Интервальную величину $\\mathbf{S}$ положим равным $K_M$ с 10%-м радиусом. \n",
    "\n",
    "[//]:  $$\\mathbf{S} = [1.2959, 1.5841] \\cdot 10^{-4}.$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Vmax =  1\n",
      "KM =    1.4400e-04\n",
      "S ⊂ [0.0001296, 0.00015841]\n"
     ]
    }
   ],
   "source": [
    "# Загрузка пакета для интервальных вычислений\n",
    "pkg load interval\n",
    "\n",
    "# Максимальная скорость реакции\n",
    "Vmax = 1\n",
    "\n",
    "# Константа Михаэлиса\n",
    "KM = 1.44e-4\n",
    "\n",
    "# Концентрация субстрата\n",
    "S = KM * infsup(0.9,1.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Для вычисления интервальной оценки $\\mathbf{v}$, выражение [(1)](#MichaelisMenten) представим двумя способами; \n",
    "в исходном виде, и с делением числителя и знаменателя на $\\mathbf{S}$:   \n",
    "\n",
    "$$\\mathbf{v}_1 = V_{max} \\ \\frac{\\mathbf{S}}{\\mathbf{S}+K_M},$$  \n",
    "\n",
    "$$\\mathbf{v}_2 = V_{max} \\ \\frac{1}{1+K_M/\\mathbf{S}}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "v1 ⊂ [0.42857, 0.57895]\n",
      "v2 ⊂ [0.47368, 0.52381]\n"
     ]
    }
   ],
   "source": [
    "# Выражение с двумя вхождениями интервальной переменной S\n",
    "v1 = Vmax * S / (S + KM)\n",
    "\n",
    "# Выражение с одним вхождением интервальной переменной S\n",
    "v2 = Vmax / (1 + KM / S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Средние величины $\\mathsf{mid}\\ \\mathbf{v}_1$ и $\\mathsf{mid}\\ \\mathbf{v}_2$ обоих выражений отличаются незначительно. В то же время радиусы результатов вычислений $\\mathsf{rad}\\ \\mathbf{v}_1$ и $\\mathsf{rad}\\ \\mathbf{v}_2$ существенно различны."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m1 =  0.50376\n",
      "m2 =  0.49875\n",
      "r1 =  0.075188\n",
      "r2 =  0.025063\n"
     ]
    }
   ],
   "source": [
    "# Середины интервальных оценок скорости реакции\n",
    "m1 = mid(v1)\n",
    "m2 = mid(v2)\n",
    "\n",
    "# Радиусы интервальных оценок скорости реакции\n",
    "r1 = rad(v1)\n",
    "r2 = rad(v2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "В первом случае результирующий интервал `v1` шире, поскольку интервальная переменная `S` входит в выражение дважды.\n",
    "\n",
    "Во втором случае, интервальная переменная `S` входит в выражение для `v2` один раз, и согласно \n",
    "основной теореме интервальной арифметики результат интервального оценивания совпадает \n",
    "с точной областью значений выражения.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Пример 2. Полное сопротивление цепи переменного тока\n",
    "\n",
    "Полное сопротивление цепи переменного тока складывается из активной и реактивной \n",
    "составляющих, общая формула \n",
    "\n",
    "<a id='Zseq'></a>\n",
    "$$Z = \\sqrt{ R^2 + (X_L - X_C)^2 }, \\hspace{3cm} (2)$$\n",
    "\n",
    "где  $R$ --- активное сопротивление, $X_L = \\omega\\cdot L$ --- индуктивное \n",
    "сопротивление, $X_C = 1/\\omega C$ --- ёмкостное сопротивление (см., к примеру, \n",
    "[[2]](#YavorskiDetlaf)). \n",
    "\n",
    "Как будет изменяться ток в цепи переменного тока с заданными \n",
    "$R$, $L$, $C$, когда круговая частота меняется в интервале $[\\,\\omega_1, \\omega_2 ]$? \n",
    "\n",
    "Для ответа на вопрос можно воспользоваться законом Ома для цепи переменного тока\n",
    "\\begin{equation*}\n",
    "I = U/Z, \n",
    "\\end{equation*}\n",
    "но до этого необходимо найти интервал изменения $Z$.\n",
    "\n",
    "Построим интервал изменения полного сопротивления $Z$ для примера цепи двумя способами: \n",
    "* сначала с помощью содержательного анализа получим точный интервал изменения $Z$;\n",
    "* затем вычислим формулу [(2)](#Zseq) в интервальной арифметике и сравним интервальную оценку с точным интервалом изменения $Z$.\n",
    "\n",
    "### Параметры электрической цепи\n",
    "\n",
    "Для расчетов конкретизируем параметры цепи.\n",
    "\n",
    "В качестве центральной частоты возьмём $f= 13.56$ МГц. Данный диапазон может \n",
    "использоваться во всем мире без лицензий,\n",
    "он имеет специальное обозначение --- [ISM](https://en.wikipedia.org/wiki/ISM_radio_band) (Industrial, Scientific, Medical). \n",
    "\n",
    "Рассмотрим практическое применение оборудования в диапазоне ISM: генератор \n",
    "плазменного источника возбуждения-ионизации пробы для элементного анализа \n",
    "с [индуктивно-связанной плазмой](https://ru.wikipedia.org/wiki/%D0%98%D0%BD%D0%B4%D1%83%D0%BA%D1%82%D0%B8%D0%B2%D0%BD%D0%BE-%D1%81%D0%B2%D1%8F%D0%B7%D0%B0%D0%BD%D0%BD%D0%B0%D1%8F_%D0%BF%D0%BB%D0%B0%D0%B7%D0%BC%D0%B0) (ИСП). \n",
    "\n",
    "Индуктор плазмотрона ИСП (диаметр и длина 20 мм, три витка) имеет индуктивность \n",
    "порядка $L=100$ нГн и реактивное сопротивление 8.52 Ом на частоте ISM. В таком случае \n",
    "ёмкость резонансного контура имеет такое же реактивное сопротивление при $C=1.38$ нФ.  \n",
    "  \n",
    "Пусть круговая частота $ \\omega$ меняется в интервале $\\ \\pm 10 \\%$ относительно \n",
    "центра  диапазона ISM. \n",
    "\\begin{equation}\n",
    "\\label{OmegaInt}\n",
    "\\mathbf{\\omega} = [\\,\\underline{\\omega}, \\overline{\\omega} ] = [7.668, 9.372] \\cdot 10^7 \\ \\text{Гц}.\n",
    "\\end{equation}\n",
    "\n",
    "Напомним, что круговая частота связана с частотой колебаний $f$ соотношением \n",
    "$\\omega = 2 \\pi f$. \n",
    "\n",
    "Активное сопротивление примем равным 1 Ом. Такой порядок имеет \n",
    "сопротивление частично ионизированного аргона в плазмотроне. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "f =  13560000\n",
      "omega ⊂ [7.6679e+07, 9.372e+07]\n",
      "L =    1.0000e-07\n",
      "C =    1.3800e-09\n",
      "R =  1\n"
     ]
    }
   ],
   "source": [
    "# Центральная частота, Гц\n",
    "f = 13.56e6\n",
    "\n",
    "# Круговая частота равна центральной ±10%, Гц\n",
    "omega = 2 * pi * f * infsup(0.9, 1.1)\n",
    "\n",
    "# Индуктивность, Гн \n",
    "L = 100e-9 \n",
    "\n",
    "# Ёмкость резонансного контура, Ф\n",
    "C = 1.38e-9  \n",
    "\n",
    "# Активное сопротивление, Ом\n",
    "R = 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Точный интервал значений $Z$ \n",
    "\n",
    "Выражение [(2)](#Zseq) для полного сопротивления $Z$ характерно тем, что индуктивное и ёмкостное сопротивления зависят от частоты противоположным образом. Это хорошо иллюстрирует график зависимости $X_L, X_C$ и $Z$ от частоты в интервале \n",
    "$\\mathbf{\\omega}.$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAJMmlDQ1BkZWZhdWx0X3JnYi5pY2MAAEiJlZVnUJNZF8fv8zzphUASQodQQ5EqJYCUEFoo0quoQOidUEVsiLgCK4qINEWQRQEXXJUia0UUC4uCAhZ0gywCyrpxFVFBWXDfGZ33HT+8/5l7z2/+c+bec8/5cAEgiINlwct7YlK6wNvJjhkYFMwE3yiMn5bC8fR0A9/VuxEArcR7ut/P+a4IEZFp/OW4uLxy+SmCdACg7GXWzEpPWeGjy0wPj//CZ1dYsFzgMt9Y4eh/eexLzr8s+pLj681dfhUKABwp+hsO/4b/c++KVDiC9NioyGymT3JUelaYIJKZttIJHpfL9BQkR8UmRH5T8P+V/B2lR2anr0RucsomQWx0TDrzfw41MjA0BF9n8cbrS48hRv9/z2dFX73kegDYcwAg+7564ZUAdO4CQPrRV09tua+UfAA67vAzBJn/eqiVDQ0IgALoQAYoAlWgCXSBETADlsAWOAAX4AF8QRDYAPggBiQCAcgCuWAHKABFYB84CKpALWgATaAVnAad4Dy4Aq6D2+AuGAaPgRBMgpdABN6BBQiCsBAZokEykBKkDulARhAbsoYcIDfIGwqCQqFoKAnKgHKhnVARVApVQXVQE/QLdA66At2EBqGH0Dg0A/0NfYQRmATTYQVYA9aH2TAHdoV94fVwNJwK58D58F64Aq6HT8Id8BX4NjwMC+GX8BwCECLCQJQRXYSNcBEPJBiJQgTIVqQQKUfqkVakG+lD7iFCZBb5gMKgaCgmShdliXJG+aH4qFTUVlQxqgp1AtWB6kXdQ42jRKjPaDJaHq2DtkDz0IHoaHQWugBdjm5Et6OvoYfRk+h3GAyGgWFhzDDOmCBMHGYzphhzGNOGuYwZxExg5rBYrAxWB2uF9cCGYdOxBdhK7EnsJewQdhL7HkfEKeGMcI64YFwSLg9XjmvGXcQN4aZwC3hxvDreAu+Bj8BvwpfgG/Dd+Dv4SfwCQYLAIlgRfAlxhB2ECkIr4RphjPCGSCSqEM2JXsRY4nZiBfEU8QZxnPiBRCVpk7ikEFIGaS/pOOky6SHpDZlM1iDbkoPJ6eS95CbyVfJT8nsxmpieGE8sQmybWLVYh9iQ2CsKnqJO4VA2UHIo5ZQzlDuUWXG8uIY4VzxMfKt4tfg58VHxOQmahKGEh0SiRLFEs8RNiWkqlqpBdaBGUPOpx6hXqRM0hKZK49L4tJ20Bto12iQdQ2fRefQ4ehH9Z/oAXSRJlTSW9JfMlqyWvCApZCAMDQaPkcAoYZxmjDA+SilIcaQipfZItUoNSc1Ly0nbSkdKF0q3SQ9Lf5RhyjjIxMvsl+mUeSKLktWW9ZLNkj0ie012Vo4uZynHlyuUOy33SB6W15b3lt8sf0y+X35OQVHBSSFFoVLhqsKsIkPRVjFOsUzxouKMEk3JWilWqUzpktILpiSTw0xgVjB7mSJleWVn5QzlOuUB5QUVloqfSp5Km8oTVYIqWzVKtUy1R1WkpqTmrpar1qL2SB2vzlaPUT+k3qc+r8HSCNDYrdGpMc2SZvFYOawW1pgmWdNGM1WzXvO+FkaLrRWvdVjrrjasbaIdo12tfUcH1jHVidU5rDO4Cr3KfFXSqvpVo7okXY5upm6L7rgeQ89NL0+vU++Vvpp+sP5+/T79zwYmBgkGDQaPDamGLoZ5ht2GfxtpG/GNqo3uryavdly9bXXX6tfGOsaRxkeMH5jQTNxNdpv0mHwyNTMVmLaazpipmYWa1ZiNsulsT3Yx+4Y52tzOfJv5efMPFqYW6RanLf6y1LWMt2y2nF7DWhO5pmHNhJWKVZhVnZXQmmkdan3UWmijbBNmU2/zzFbVNsK20XaKo8WJ45zkvLIzsBPYtdvNcy24W7iX7RF7J/tC+wEHqoOfQ5XDU0cVx2jHFkeRk4nTZqfLzmhnV+f9zqM8BR6f18QTuZi5bHHpdSW5+rhWuT5z03YTuHW7w+4u7gfcx9aqr01a2+kBPHgeBzyeeLI8Uz1/9cJ4eXpVez33NvTO9e7zofls9Gn2eedr51vi+9hP0y/Dr8ef4h/i3+Q/H2AfUBogDNQP3BJ4O0g2KDaoKxgb7B/cGDy3zmHdwXWTISYhBSEj61nrs9ff3CC7IWHDhY2UjWEbz4SiQwNCm0MXwzzC6sPmwnnhNeEiPpd/iP8ywjaiLGIm0iqyNHIqyiqqNGo62ir6QPRMjE1MecxsLDe2KvZ1nHNcbdx8vEf88filhICEtkRcYmjiuSRqUnxSb7JicnbyYIpOSkGKMNUi9WCqSOAqaEyD0tandaXTlz/F/gzNjF0Z45nWmdWZ77P8s85kS2QnZfdv0t60Z9NUjmPOT5tRm/mbe3KVc3fkjm/hbKnbCm0N39qzTXVb/rbJ7U7bT+wg7Ijf8VueQV5p3tudATu78xXyt+dP7HLa1VIgViAoGN1tubv2B9QPsT8M7Fm9p3LP58KIwltFBkXlRYvF/OJbPxr+WPHj0t6ovQMlpiVH9mH2Je0b2W+z/0SpRGlO6cQB9wMdZcyywrK3BzcevFluXF57iHAo45Cwwq2iq1Ktcl/lYlVM1XC1XXVbjXzNnpr5wxGHh47YHmmtVagtqv14NPbogzqnuo56jfryY5hjmceeN/g39P3E/qmpUbaxqPHT8aTjwhPeJ3qbzJqamuWbS1rgloyWmZMhJ+/+bP9zV6tua10bo63oFDiVcerFL6G/jJx2Pd1zhn2m9az62Zp2WnthB9SxqUPUGdMp7ArqGjzncq6n27K7/Ve9X4+fVz5ffUHyQslFwsX8i0uXci7NXU65PHsl+spEz8aex1cDr97v9eoduOZ67cZ1x+tX+zh9l25Y3Th/0+LmuVvsW523TW939Jv0t/9m8lv7gOlAxx2zO113ze92D64ZvDhkM3Tlnv296/d5928Prx0eHPEbeTAaMip8EPFg+mHCw9ePMh8tPN4+hh4rfCL+pPyp/NP637V+bxOaCi+M24/3P/N59niCP/Hyj7Q/Fifzn5Ofl08pTTVNG02fn3Gcufti3YvJlykvF2YL/pT4s+aV5quzf9n+1S8KFE2+Frxe+rv4jcyb42+N3/bMec49fZf4bmG+8L3M+xMf2B/6PgZ8nFrIWsQuVnzS+tT92fXz2FLi0tI/QiyQvpTNDAsAAAAJcEhZcwAACxMAAAsTAQCanBgAAAAddEVYdFNvZnR3YXJlAEdQTCBHaG9zdHNjcmlwdCA5LjI2WJButwAAGMBJREFUeJzt3T9MI/mhwPEhSg1LyqedfdFGWhemIykwXY6TcBXFEst1AenwlaRYxBXhitsmHDSUByctxSsWVkFKE1Y68ro1Uu6ozi5IQYFPr4n0WCOl5hXW83mN+eux+dnz+VTGOzfM/m52vszPw8zQxcVFBAAP7WcPvQEAEEWCBEAgBAmAIAgSAEFIIEi1Wq1arTa+LJfLtVqt89UCkCoJBGl3d/f169f113Nzczs7O/Pz86VSqfM1A5AenQapUCisr6/XX+/v78dx/PLly42Njc3NzY63DYAU6TRIe3t7n376af11uVweGxuLoiiO48PDw043DYA0SfiihjiO6y8mJiaSXTMAg+3nya6ucXXDVWdImUwm2e8IQMPx8fFDb8L9JRmksbGxcrkcRVG1Ws1ms1ct1tfjFaahoSG3gOoGA9slBrZL+v0n/iSDlM/nt7a21tbWDg8PFxYWElwzAAMv+Z9TSqVSHMeND5NaZDIZZ0iJ8/NmlxjYLjGwXdLvB9iEP0OKoiiXyyW+TgAGnlsHDQI/bHaJge0SA0tbggRAEAQJgCAIEgBBECQAgiBIAARBkAAIgiABEARBAiAIggRAEAQJgCAIEgBBECQAgiBIAARBkAAIwuAHaejFfz/0JgBws+Qf0Beai/XftjTpYv23D7UxAFxl8IMUXSqQPgEEKBVBaqFPwAP658x/RFH07M3/PPSGBCeNQWqhT0D31PPTTIquIkit9Am4N/nphCDd4Po+XV4ASA/5SZYg3c3l/DiFgpSQn24TpE6Z4oOBJD+9J0gJ0yfoR/ITAkHqLn2CAMlPmASpp1wiAb0nP/1CkB6SSyQgWZfbE8lP/xCksJjig9tz6jNgBClo+gQN8jPwBKmf+AiK9JCfFBKkPuYjKAaG/BAJ0oBxCkVfkB/aEqRB5hSKEMgPtyRI6eIqCbpNfrg3QUo1U3x0SH5IkCDxkxun+NouQ3rID10lSFzHp1Cp5ZYH9J4gcTdm+QaSUx9CIEh0xCxfP5IfwiRIJMwsX2jkh34hSHSdWb5ekh/6lyDRa2b5EiQ/DBJB4uFJ1C3JD4NNkAiRD6Ii+SF9BIn+MPAfRMkPCBJ9qd9n+eQHLhMkBkTIiZIfuA1BYmA9VKLkB+5HkEiRbiRKfiApgkSq3TVR8gPdI0jwgUZ+Gu3558xPf5r5z/+Kwr5cAvqXIEEU3frU5yKKopAul4BBIkikUYczbyFf0Qf9S5AYfD344Oc2iWq7GNAgSAyaQK47aNseJ1JwDUGivwWSn1sy1wfXECT6SX/l5zbM9UGDIBGuwcvPbZjrI7UEiVCkMz+3ZK6PNBAkHob8dMhcH4NHkOgF+ekBc330O0EiYZfbE8nPw3EiRR8RJDri1KfvOJEiWILEHcjPoHIiRQgEiSvJT5rd8kTqqiXhHpIPUrlcjuN4ZGQk8TXTVfLDjUz30VVJBqlWq83Pz2ez2VKpVCwWZ2dnE1w5yZIfkmK6j6QkGaS3b99OTEwsLS3VarXFxUVBCof80Eum+7ifJIOUzWZ3dnZKpVKlUslmswmumTuRHwKkUtwoySDFcTw8PLy/v1+pVK45PRoaGmq8vri4SHAD0kl+6F8q1bnmI+qzZ88ecEs6l2SQNjc38/l8PUUfffTRVU0SoU7IDwNPpe6k+YiayWQecEs657LvoMkP1LnALw2SDNInn3wyNzd3enpaqVSmp6cTXHNKyA/ciQv8BsxQ4hNopVIpjuM4jtv+aSaTOT4+TvY79in5gd5oW6loEEPV7wfY5Kfscrlc4uvsd+43Cg/oqvA4nQqNz5CS59QH+oKrJ0IjSJ2SHxgkKvWABOlu5AdSSKV6Q5CuIz/AVW5fqasWpoUg/UR+gA7d/gIKibosvUGSH6Bn5Oc20hIk+QEI3OAHqZ4i+QEI3OAHSYoA+sLPHnoDACCKBAmAQAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQBAmAIAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQBAmAIAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQBAmAIAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQBAmAIAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQBAmAIAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQBAmAIAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQkg9StVqtVquJrxaAwfbzZFe3srISRVG1Ws3n87Ozs8muHIABluQZUqlUiqLo5cuXGxsb9dcAcEtJniG9e/fuyZMnOzs7URRtbGwkuGYABl7CnyG9fv06iqLT09O5ubmrlhlqkux3B0ibQTqitjlDqtVqb9++PT09bbyztLR0y9VNT0/XPzoqFApXLXNxcXHHjQSgveYjaiaTecAt6VybM6TFxcX7fQL05MmTxutarXb/jQIgfdqcIVUqle++++4e65qdnS0UCmtra5VKpVgsdrxtAKTI0OUJtLW1tUePHi0sLNxvjaVSKY7jOI7b/mkmkzk+Pr7fmgG4Rr8fYNtfZbe+vr6+vt748k5/w1wu1+lGAZA+bYL0zTff/OMf/xgZGen91gCQWm0uapiYmOj9dgAwwE5OTk5OThpfHh0dnZ2dtSzT5gwpjuOpqalsNtt4Z3t7uztbCEAqnJ2dzczMHB0djY6OHhwcLC8vHx0dtSzTJkizs7P5fL4nWwhAKoyPj3/22WfFYnFzc3N5eXlzc/PyMm2m7K66QA4A7m15ebl+nvT8+fPx8fHLC7Q5Q/riiy+iKHr8+HHjHRfOAdC5jz/++PPPP297ehS1DdLbt2/7+kp2AAJ0cnLy9ddff/3118Vi8eDgoP7m5ubm999/H0XRr3/96zZBmp6eLpfLY2NjPd1SAB7C0Iv/TnydF+u/vfxmsVhcXV2dmZn59ttvV1dXl5eXDw4Ovv/++/oJ08zMzAd3aqjfovv8/LxarXbpKrt+/0VigGCFfIBdXV39/vvv37x5E0XR2dnZ+Pj4wcHB119//fHHH09NTdWX+eAMqXEDun/961///ve/f/nLX/Z2gwEYTMvLy43Xo6Ojzb+T1PBBkHK53Nra2u7ubhRFw8PD5+fnz58/v/2zJwDg9n71q18dHR3Vz5CWl5c/mLJbW1s7PDzc2NioX/ldrVbX19cfP36cYJNCPqME6Gv9eICdmZl5+vRpFEUnJycfBOmjjz7a29trvotdrVYrFAotb3aiH8cLoC/06QG2Pn339OnTn34xtlwux3HcEp6RkZE4jiuVSq83EIB0ePr0af0k6acgxXFcrVYvL3p+fj48PNy7TQMglX4KUv1kaGdnp/mP61/6nSQAuu2Dq+w2Njbm5+f39/frv4T0448/lsvljY2NB9o2AFLkg5urjoyM7O3tzc7O1r/M5XJ7e3tOjwDogTa3Dsrn8x4/AUCPtXn8BAD0niABEARBAiAIggRAEG4IUqlUymQyW1tbLb+fBADJ+iBIa2trtVqt+Z1cLvfpp58uLCw0rgUHgG74IEjffPNNoVAol8stb/Z2kwBIow+CNDExUSwWFxcXTdABkKCTk5Pmh/IdHR2dnZ21LNP6GdLs7OzGxsbOzs7KykrXNxCAdDg7O5uamqpH6ODgoFgsjo6OtizT5qKGsbGxV69enZ+fFwqFtvf/BoA7GR8f/+yzz4rF4tnZ2fLy8ubm5uVl2l9lNzIysrGxkc/nC4VClzcSgFRYXl4+OzubmZl5/vz5+Pj45QU+CFLLLewWFhY2NjYmJia6u40ApMPHH3/897//fWZmpu2ffvAI8x7o0yfsAoQv8APsycnJ1NTU559/vru7e3BwEEXR0dFRFEWNs6U2d/sGID3+OfMfia/z2Zv/ufxmsVhcXV2dmZn59ttvV1dX6zN4zQsIEkCqtY1H4lZXV0dHR+uTdZubm+Pj45cn7gQJgK5bXl5uvB4dHa3/TlLzbyZF7mUHQCBuOEPK5XIhf0QGQF8rFotPnz6Nomh1ddWUHQAPY2pqqnnW7ropu7W1tfpMXalUunzTVQBI0JVB2tnZOT8/rz91IpfLbWxsfPHFFz3cMADS5cognZ6eNt+4IY7jiYmJUqnUk60CIHWuDNKjR49a7qxaqVTiOO7+JgGQRlde1PD8+fPFxcXT09PJyclarfb27dtsNitIAHTJDfey29nZOT09jaJocnIyl8t1/v0Cv9USQP/q9wPsB1N2i4uLtVqt+Z1cLnd4eLi0tJRIjQDgKh8E6fHjx4VCoXHlwtbW1tzcXMszKQCgGz74DGlpaWlycnJlZWV6erpSqURRtL297XMjAHqg9aKGXC734sWLP/7xj8PDw69evVIjAHrjgyDVarUvvviiWq3+5S9/qVar8/Pzz58/X1paeqiNAyA9PvgMaWpq6vHjx3t7e2NjY/l8/uDgoFKpFAqFh9o4ANLjgzOkjY2N5qvpRkZGtre3t7a2er5VAKTODb+HlLh+v0weIFj9foD1+AkAuu7k5KT5SROjo6Pj4+MtywgSAF13cnLy7bff1l8fHR2Njo6+efOmZRlBAqDrpqampqamoig6Ojp68+bNwcHB5WWue0BfFEWlUimTyWxtbdWf1AcA93Z2dlYsFldXV+uPLW9xwxlSLpfr64/IAAjH8vLy1NTUzMxM2z+9w5Td2tqaX5IFGDD/+9c2Jysd+sXvTi6/ubm5eXJy0nayru4OQarf3Q6AQdI2Hok7Ojr685//fHR0dM0yLmoAoOuWl5ejKGpM1o2Pj6+urrYsc2WQSqXSu3fvmt9peaI5ANzSNTN1DVcGKY7jycnJ5ndavgSABF35xNg4jnO5XBzH6+vruf/3EFsIQCp4YiwAQfDEWACC0JUnxlar1eHh4ZGRkSS2EIBUSP6JsbVarVAofPfdd4luJwADLvknxq6vrw8PDzcujgCA20j4ibFbW1tPnjypVqvm6wC4kw/OkNpe2L2wsHDLdZXL5XK5fOPyQ01uuWYA2hqkI+oNj5+4k62trcePH6+trVWr1ZWVlavu7HDRJMHvDpBCg3RETfJedgsLC+fn51EUHR4e5vP54eHhBFcOwGBLMkhjY2P1F8PDw27rAMCdJDll17C9vd2N1QIwwLoSJAC4K0ECIAiCBEAQBAmAIAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQBAmAIAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQBAmAIAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQBAmAIAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQBAmAIAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQBAmAIAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQBAmAIAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQBAmAIAgSAEEQJACCIEgABEGQAAiCIAEQBEECIAiCBEAQkg9SuVyuVquJrxaAwfbzBNdVq9Xm5+ez2Wy1Ws1ms0tLSwmuHIDBluQZ0u7u7sTExMuXL7e3t9++fZvgmjvxv399+tCbAMDNkjxDmp6err+o1WoJrrZDv/jdSXOTfvG7kwfcGACukmSQ4jiOoqhUKq2vrxeLxasWGxoaary+uLhIcAOu0hyhRpyUCRgAzUfUZ8+ePeCWdC7JIEVRtLa29uOPP25sbNTj1FZvInSVRoeUCRgAzUfUTCbzgFvSuSSDtLOzU69RguvsHmUCCEqSQapf8D03N1f/cnt7O8GVd48yAT3jOHONoR5PoGUymePj415+x3twEQSQoJ4dUvriAHuNhD9DGgwuggA65OfaexCkG5jQA25JhDokSLelTMBlIpQgQbqzy2WK7IiQJv7td4kg3Z+PmiA9RKgHBCkZJvRg8IhQjwlSwkzoQV/zL/cBCVK3mNCDfiFCgRCkXnDaBKHxjzFAgtRTbU+bIv8eoCf8owucID0Yc3rQAyLURwQpCOb0ICktD4n276iPCFJYzOnBPfjHMhgEKVzm9OAaIjR4BKk/mNMDc3EDT5D6jDk9UsVOniqC1MfEiYFkZ04tQRoQ4kT/MhdHnSANIHEicApEW4I04MSJQNj9uJEgpchVcYocIOgC+xh3JUgp1XJ08NMrnVMgOiRIRJGTJ+7FrkKyBIlW15w8Xf5TUsXOQFcJEjfQp9Rq+X8d+d9NlwkSd6NPg0p+eHCCREf0qU/JDwESJJJ0fZ8uL0Bv+B9BXxAkuujyUe/ykbHtYtybEaZ/CRI91fbI6Bh6P8aNASNIPLzbV+qa5QeVcSA9BIlAXXO0veYYff1/GKDr/y5Rv/11oBOCRP+5/hh94yH+Huu8jYf6vjAwBIlBc79D/P1y0vn3BRoECaJITiAAP3voDQCAKBIkAAIhSAAEQZAACIIgARAEQQIgCIIEQBAECYAgCBIAQRAkAIIgSAAEQZAACIIgARAEQQIgCIIEQBAECYAgCBIAQRAkAIIgSAAEQZAACIIgARAEQQIgCIIEQBAECYAgCBIAQRAkAIIgSAAEQZAACIIgARAEQQIgCIIEQBCSD1K5XK7VaomvlmsMDQ099CYMJgPbJQaWthIO0tzc3M7Ozvz8fKlUSnbNAAy2nye4rv39/TiOX758Wa1WV1ZWcrlcgisHYLAleYZULpfHxsaiKIrj+PDwMME1AzDwEp6yi+O4/mJiYiLZNQMw2JKcsouiqFqt1l9cc4aUyWSS/aY8e/bMqHaDge0SA0tbSQZpbGysXC5HUVStVrPZbNtljo+PE/yOAAyMoYuLiwRXVygUJiYmDg8PFxYW8vl8gmsGYLAlHKQoikqlUhzHjQ+TAOA2kg8SANxDQLcOcouHq9xjZAzmVRIcGYN8FXtsglK1x97/ooadnZ39/f3660ql8urVq/ovIUVR9Jvf/KZ+UUM2m11aWrrN2ubm5uI4rlQqL168yOVy16w8VWq12vz8fDabLZVKxWJxdnb2Nv9V82BWq1UjWXfNYNpjk1Kr1RYXF7PZ7J0+SLbHtnXNYHa4x4Y7yBcd++GHH7766qvGl6enp3/4wx9alnn//v27d+9++OGHtmv429/+9qc//antf9uy8rR5/fp1/a///v37+shcP5IXVw9mykfyot1g1tljE/TVV1+9fv364uLi/fv3v//97y/ssR24PJh1l/e6gRnkBKbs1tfXi8Vi48tqtToyMrKysrK2tlY/PazVaoVC4d27d1tbW2tra/XFdnZ26teIR9fe4qFl5WlT/+GoVCrt7u5ms9m2IxndbjBTPpLRpcFsvG+PTdbw8HAURSMjI5VKxR7boebBbLzZsscO0iB3+ntIOzs7k5OTIyMjjXdqtdrjx48nJycrlcri4uL29vbu7m5jhuSjjz6qn2Cenp42X4nX9hYPl1eeNnEcDw8P7+/vVyqV2dnZtiMZ3WIwjWR0aTAb79tjEzQ5Obm+vh5F0du3b4eHh+2xnWgZzMb7LXvs5OTkwAxyAkF69epV8zv5fL4+15nL5erTlO/fv9/f36+/juO4PtterVYPDw+Hh4dfvHgRXXGLh8srT5vNzc18Pt/Y1aanp5tHMvr/T/JuHEwjGV0azEaT7LEJyuVyX3755f7+/uzsbKlUahnJyB57Fy2D2Xi/ZY8dpEHuKEjlcjmbzbYEdmtrK5vNNt/q+9GjR7Ozs/V//1tbW/XXa2trk5OT9cWq1erlWzy0XXnKtYxkFEW3GUwjeT17bIL29/fPz8+Xlpaq1Woul7PHdqJlMBvvt+yxgzTIHQVpf3+/cW1GqVSan58/Pj6emJhYXFycnp5uTIw8f/58cXHx9PS0Uqm0veomn8/XZz/rF5NcXnlqffLJJ3Nzc/Whm56evnEko3aDaSTrWgbTHtsNuVyuUCicnp4eHh5++eWXcRzbY++tZTCv2mOnp6cHZpC79Yuxl+/XcOMdHNzi4SotI3ObgTKYV7lqZOyxSanVapVKpfmnb3vsvV0ezIaBPCy4UwMAQQjoTg0ApJkgARAEQQIgCIIEQBAECYAgCBIAQRAkuEEmk7nmSyApggRAEDq9uSr0qXK5vL+//8knn0RR9Pr163w+f6c7qTQ/ke/2D0kDriFIpFH96bF7e3v1ByDl8/n5+fmDg4Or7jXZ/JiZusbtLFdWVrq9tZASpuxIo/oz+uI4rt+3f2xsLJvN7u7uXrX8ZJPm9+s3V3Z6BIlwhkQavX///pZv1jXf/L+hVCrt7++H9kQZ6F/OkEijJ0+e1F+MjY01Pjqqv9n8JLRrlMvllZWVjY2N0J4oA/3L3b5JqUKhMDEx8ejRoyiK3r9/f3h4uLe3F0VRJpM5Pj5uXrLlnfqXi4uL5XK5cSf/7e3t3m06DChBIr2q1Wq9KI0XURRtbW01nrkH9JIpO9KrEaHmR5YNDw8/0OZA2jlDAiAI/wcFHDIlxzXY5gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Сетка значений в интервале omega\n",
    "omega_vals = omega.inf : rad(omega)/20 : omega.sup;\n",
    "\n",
    "# Значения ёмкостного сопротивления\n",
    "XC = 1 ./ (omega_vals .* C);\n",
    "\n",
    "# Значения индуктивного сопротивления\n",
    "XL = omega_vals .* L;\n",
    "\n",
    "# Значения полного сопротивления\n",
    "Z = (R^2 + (XL - XC) .^ 2) .^ 0.5;\n",
    "\n",
    "# График\n",
    "plot(omega_vals, XC)\n",
    "hold on\n",
    "plot(omega_vals, XL)\n",
    "plot(omega_vals, Z)\n",
    "xlabel('\\omega, Hz')\n",
    "ylabel('X_L, X_C, Z, Ohm')\n",
    "legend({'X_C','X_L','Z'},\"location\", \"eastoutside\")\n",
    "legend boxoff"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Определим границы изменения величины $Z$. \n",
    "\n",
    "В точке резонанса  $X_L = X_C$ суммарное \n",
    "сопротивление $LC$-части схемы равно нулю, а полное сопротивление равно активному, \n",
    "т.е. $Z = R$. Таким образом, \n",
    "$$\\min_{\\omega} Z = R = 1.$$\n",
    "\n",
    "Величины индуктивного $X_L$ и ёмкостного $X_C$ сопротивлений зависят от частоты \n",
    "противоположным образом и монотонны. Поэтому максимальные значения величины \n",
    "$\\Delta = | X_L - X_C | $ достигаются на одном из краёв диапазона $\\mathbf{\\omega}$: \n",
    "$$\n",
    "\\Delta_{max} = \\max \\left\\lbrace \\Delta( \\underline {\\mathbf{\\omega}}), \\Delta( \\overline {\\mathbf{\\omega}}) \n",
    "\\right\\rbrace =  \\max \\left\\lbrace 1.7822, 1.6401  \\right\\rbrace = 1.7822. $$\n",
    "\n",
    "Соответственно, \n",
    "$$\\max_{\\omega} Z = \\sqrt{R^2 + \\Delta_{max}^2}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Zmin =  1\n",
      "DeltaLeft =  1.7822\n",
      "DeltaRight =  1.6401\n",
      "DeltaMax =  1.7822\n",
      "Zmax =  2.0435\n",
      "Zexact ⊂ [1, 2.0436]\n",
      "rexact =  0.52177\n"
     ]
    }
   ],
   "source": [
    "# Нижняя граница полного сопротивления\n",
    "Zmin = R\n",
    "\n",
    "# Delta для левой границы интервала omega\n",
    "DeltaLeft = abs(XL(1) - XC(1))\n",
    "\n",
    "# Delta для правой границы интервала omega\n",
    "DeltaRight = abs(XL(end) - XC(end))\n",
    "\n",
    "# Максимум Delta - максимум значений на границах интервала omega\n",
    "DeltaMax = max(DeltaLeft, DeltaRight)\n",
    "\n",
    "# Верхняя граница полного сопротивления\n",
    "Zmax = (R^2 + DeltaMax^2) ^ 0.5\n",
    "\n",
    "# Точный интервал изменения полного сопротивления\n",
    "Zexact = infsup(Zmin, Zmax)\n",
    "\n",
    "rexact =  rad(Zexact)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Окончательно имеем точный интервал \n",
    "$$\\mathbf{Z} = [1, 2.044].$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Вычисление интервального выражения\n",
    "\n",
    "Вычислим также интервал $\\mathbf{Z}$ по формуле [(2)](#Zseq), представленной двумя способами:\n",
    "\n",
    "$$\\mathbf{Z}_1 = \\sqrt{R^2 + \\left( \\mathbf{\\omega} L - \\frac{1}{\\mathbf{\\omega} C} \\right) ^2},$$ \n",
    "\n",
    "$$\\mathbf{Z}_2 = \\sqrt{R^2 + \\mathbf{\\omega}^2\\left( L - \\frac{1}{\\mathbf{\\omega}^2 C} \\right)^2}.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Z1 ⊂ [1, 2.0436]\n",
      "Z2 ⊂ [1, 2.3968]\n",
      "r1 =  0.52177\n",
      "r2 =  0.69838\n"
     ]
    }
   ],
   "source": [
    "# Выражение с двойным вхождением omega\n",
    "Z1 = (R^2 + (omega * L - 1 / (omega * C))^2)^0.5\n",
    "\n",
    "# Выражение c двойным вхождением квадрата omega\n",
    "Z2 = (R^2 + omega^2 * (L - 1 / (omega^2 * C))^2)^0.5\n",
    "\n",
    "# Радиусы интервалов\n",
    "r1 = rad(Z1)\n",
    "r2 = rad(Z2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Нижние величины интервалов $Z_1, \\ Z_2$ соответствуют частотному резонансу $X_L=X_C$ в середине диапазона $\\mathbf{\\omega}$, верхние относятся к краям выбранного диапазона, иначе --- расстройке резонанса. \n",
    "\n",
    "Значение выражения $Z_1$ в точности совпадает с точным интервалом изменения значений полного сопротивления, полученным аналитически. Радиус значения выражения $Z_2$ больше в связи с дополнительным умножением на величину $\\omega^2$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a id='refs'></a>\n",
    "## Публикации\n",
    "\n",
    "1. <a id='SharyBook'></a> Шарый С.П. *Конечномерный интервальный анализ.* -- ФИЦ ИВТ: \n",
    "Новосибирск, 2021.  \n",
    "http://www.nsc.ru/interval/Library/InteBooks/SharyBook.pdf \n",
    "\n",
    "2. <a id='YavorskiDetlaf'></a> Яворский Б.М., Детлаф А.А. *Справочник по физике для инженеров и студентов ВУЗов*. -- Москва: Наука, 1968 (*а также последующие издания*). \n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Octave",
   "language": "octave",
   "name": "octave"
  },
  "language_info": {
   "file_extension": ".m",
   "help_links": [
    {
     "text": "GNU Octave",
     "url": "https://www.gnu.org/software/octave/support.html"
    },
    {
     "text": "Octave Kernel",
     "url": "https://github.com/Calysto/octave_kernel"
    },
    {
     "text": "MetaKernel Magics",
     "url": "https://metakernel.readthedocs.io/en/latest/source/README.html"
    }
   ],
   "mimetype": "text/x-octave",
   "name": "octave",
   "version": "4.2.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}