\documentclass[11pt]{article}
%\usepackage[a5paper]{geometry}
\usepackage[left=20mm, right=20mm, top=20mm, bottom=20mm, includefoot]{geometry}
%\geometry{left=20mm, right=20mm, top=20mm, bottom=20mm, includefoot}

\usepackage[cp1251]{inputenc}

\usepackage[russian]{babel}
\usepackage[intlimits]{amsmath}
\interdisplaylinepenalty=2500
\usepackage{amsfonts}
\usepackage{amsthm}
\usepackage{amssymb}
\usepackage{euscript}
\usepackage{amsmath}
\usepackage{latexsym}
\usepackage{graphicx}

\DeclareMathOperator{\sech}{sech}
 %нужен, чтобы вставлять рисунки
%\usepackage{caption} %нужен, чтобы убрать автоматическую нумерацию рисунков
\newcommand{\RNumb}[1]{\uppercase\expandafter{\romannumeral #1\relax}}

\begin{document}

\begin{titlepage}
	\begin{center}

			Московский государственный университет им. М.В. Ломоносова
		\vspace{0.25cm}

		Механико-математический факультет

		Кафедра теоретической механики и мехатроники
		\vfill



		\textsc{Отчет}\\[5mm]

		{\LARGE По задаче 5 вычислительного практикума №2}
		\bigskip

		Жучков Алексей
		\vfill
		4 курс, группа 422
	\end{center}
	\vfill



	\begin{center}
		Москва, 2017 г.
	\end{center}
\end{titlepage}
\begin{center}

\end{center}
\section{Постановка задачи.}
Рассмотрим задачу Лагранжа с фиксировнным временным отрезком и без ограничений вида "меньше или равно"
$$
\int\limits_0^{\pi/2} u^2 dt + x^2(0) \to \inf, \qquad \ddot x + x \exp(-\alpha x) = u,$$

$$
\dot x(0)=x(\pi/2)=0, \dot x(\pi/2)=1, \qquad \alpha \in [0.0, 25.0].
$$

\section{Формализация задачи как задачи оптимального уравнения.}

Введем новую переменную $y=\dot x$.  Тогда система перепишется в следующем виде:

\begin{equation}\label{lag}\left\{
\begin{aligned}
&\dot x = y,\\
&\dot y = u - x \exp(-\alpha x),\\
&y(0) = 0,\\
&y(\pi/2) = 1,\\
&x(\pi/2)=0,\\
&B_0 = \int_0^{\pi/2} u^2 dt +x^2(0) \to \inf,\\
&\alpha \in [0.0, 25.0]
\end{aligned}\right. \end{equation}

\section{Системы необходимых условий оптимальности.}

Запишем функции Лагранжа и Понтрягина

$$
\mathcal L = \int \limits_0^{\pi/2} Ldt + l,
$$

\begin{tabbing}
 лагранжиан \quad \= $L = \lambda_0 u^2 +p_x (\dot x - y) + p_y(\dot y - u + x \exp(-\alpha x))$\\
терминант \quad \> $l = \lambda_0 x^2(0)+ \lambda_1 y(0) + \lambda_2 (y(\pi/2)-1) + \lambda_3 x(\pi/2)$\\
\> $H = p_x y + p_y(u - x \exp(-\alpha x)) - \lambda_0 u^2$\\
\end{tabbing}

Применим к задаче принцип максимума Понтрягина. Выпишем необходимые условия оптимальности.

\begin{enumerate}
\item Уравнения Эйлера-Лагранжа:
\begin{tabbing}
$\dot p_x = - \frac{\partial H}{\partial x}$\quad \= $ \dot p_x=  p_y \exp(-\alpha x)-\alpha p_y x \exp(\alpha x)$ \\
$\dot p_y = - \frac{\partial H}{\partial y}$ \quad \>   $\dot p_y= -p_x$
\end{tabbing}
\item Условие оптимальности по управлению: $u = \arg abs \max \limits_{u\in \mathbf R} H(u)$:
$u = \arg abs \max \limits_{u\in \mathbf R}(p_y u - \lambda_0 u^2) = \frac{p_y}{2\lambda_0}$. В самом деле, $H$ как функция $u$ - парабола с ветвями, направленными вниз. Максимум она достигает в точке вершины при $u=\frac{p_y}{2\lambda_0}$.
Считаем $\lambda_0 \neq 0$.
\item Условия трансверсальности по $x, y$,  $p_x(t_k)=(-1)^k \frac{\partial l}{\partial x(t_k)}$, $p_y(t_k)=(-1)^k \frac{\partial l}{\partial y(t_k)}$, где $k=0,1$, $t_0 = 0$, $t_1 = \pi/2$:
    $$\left\{
    \begin{aligned}
    p_x(0) = \frac{\partial l}{\partial x(0)} &= 2 \lambda_0 x(0),\\
    p_y(0) = \frac{\partial l}{\partial y(0)} &= \lambda_1,\\
    p_x(\pi/2) = -\frac{\partial l}{\partial x(\pi/2)}&= - \lambda_3,\\
    p_y(\pi/2) = -\frac{\partial l}{\partial y(\pi/2)}& = -\lambda_2.\\
    \end{aligned}
\right.
$$
\item Условия стационарности по $t_k$ - в данной задаче отсутствуют, так временной промежуток фиксирован.
\item Условия дополняющей нежёсткости - в данной задаче отсутствуют, так как нет условий вида "меньше либо равно".
\item Условие неотрицательности: $\lambda_i \ge 0$.
\item Условие нормировки (множители Лагранжа выбираются с точностью до положительной константы).
\item НЕРОН (все множители Лагранжа одновременно не равны нулю).

\end{enumerate}
\section{Анормальный случай}
 Рассмотрим анормальный случай $\lambda_0 = 0$. Тогда если на отрезке $[0, \pi/2]$ существует интервал $(t_1,t_2)$ такой, что при $t\in (t_1, t_2)$ имеем $p_y(t) \ne 0$, то из условия оптимальности управления получим $u\to \infty$, и такой управляемый процесс не принадлежит допустимому классу процессов. В таком случае имеем $p_y(t)\equiv 0$ на всем отрезке $[0, \pi/2]$. Тогда из условий трансверсальности имеем $\lambda_1=\lambda_3=0$. Кроме того, из уравнений Эйлера-Лагранжа получим, что $\dot p_x = 0$, то есть $p_x(t)\equiv c$, где $c$ - константа, то есть из условий трансверсальности: $p_x(0) = 2 \lambda_0 x(0) = 0 = p_x(\pi/2) = -\lambda_2$. В результате получили, что все множители Лагранжа нулевые, что противоречит условию НЕРОН.

Итак, доказано, что $\lambda_0 \ne 0$. В силу однородности функции Лагранжа по множителям Лагранжа выбираем следующее условие нормировки: $$\lambda_0 = \frac{1}{2}$$.
Тогда из условия 2 определяется управление и начальное условие для сопряженной переменной $p_x$: $$u=p_y, \qquad p_x(0)= x(0)$$ .
\section{Краевая задача.}
Итак, на основе принципа максимума Понтрягина мы свели задачу к следующей краевой задаче:

\begin{equation}\label{system}
\left\{
\begin{aligned}
&\dot x = y,\\
&\dot y = p_y - x \exp(-\alpha x),\\
&\dot p_x = p_y \exp(-\alpha x)-\alpha x p_y \exp(-\alpha x),\\
&\dot p_y = -p_x,\\
&y(0) = 0,\\
&y(\pi/2) = 1,\\
&x(\pi/2)=0,\\
&p_x(0)=x(0),\\
&\alpha \in [0.0, 25.0]
\end{aligned}\right.\end{equation}
\section{Численное решение краевой задачи методом стельбы.}
Решаем краевую задачу численно методом стрельбы. В качестве параметром пристрелки выбираются недостающие начальные условия для постановки задачи Коши, а именно $x(0)=\alpha_1, p_y(0) = alpha_2$.
Выбрав их каким-либо образом, решаем задачу Коши на отрезке $[0,\pi/2]$, получаем четыре функции $x(.)[\alpha_1,\alpha_2],\quad y(.)[\alpha_1,\alpha_2],\quad p_x(.)[\alpha_1,\alpha_2],\quad p_y(.)[\alpha_1,\alpha_2]$ зависящих от выбора конкретных значений $\alpha_1, \alpha_2$. Решаем задачу Коши, состоящую из системы дифференциальных уравнений и начальных условий \ref{system}. с добавлением начальных условий $x(0)=\alpha_1, p_y(0) = \alpha_2$.
 численно явным методом Рунге-Кутты, основанным на расчетных формулах Дормана-Принса 5(4) DOPRI5 (см. [3])с автоматическим выбором шага .

Итак, необходимо подобрать $\alpha_1, \alpha_2$ так, чтобы выполнялось следующее условие:
$$\begin{aligned}
x(\pi/2)[\alpha_1,\alpha_2]&=0\\
y(\pi/2)[\alpha_1,\alpha_2]-1& =0
\end{aligned}$$

Соответственно, в качестве вектор-функции невязок берем функцию $$X(\bar \alpha) = \begin{pmatrix} X_1\\X_2\end{pmatrix}=\begin{pmatrix} x(\pi/2)[\alpha_1,\alpha_2]\\y(\pi/2)[\alpha_1,\alpha_2]-1\end{pmatrix}$$
Таким образом, в результате выбора вычислительной схемы метода стрельбы, решение краевой задачи свелось к решению системы двух алгебраических уравнений от двух неизвестных. Корень $\bar \alpha$ системы алгебраических уравнений $X(\bar \alpha)=0$ находится методом Ньютона с модификацией Исаева-Сонина (см. [2], [4]) по формулам:
\begin{equation}\begin{split}
X'(\bar\alpha_n)\bar\alpha_{n+1}&=\bar\alpha_n - \gamma  X(\bar\alpha_n)\\
X'(\bar\alpha_n)&=\begin{pmatrix} a & b \\ c& d \end{pmatrix}\\
a&=\frac{X_1(\alpha_{n_1} +\varepsilon,\alpha_{n_2})- X_1(\alpha_{n_1},\alpha_{n_2})}{\varepsilon}\\
b&=\frac{X_1(\alpha_{n_1},\alpha_{n_2}+\varepsilon)-X_1(\alpha_{n_1},\alpha_{n_2})}{\varepsilon}\\
c&=\frac{X_2(\alpha_{n_1} +\varepsilon,\alpha_{n_2}- X_2(\alpha_{n_1},\alpha_{n_2})}{\varepsilon}\\
d&=\frac{X_2(\alpha_{n_1},\alpha_{n_2}+\varepsilon)-X_2(\alpha_{n_1},\alpha_{n_2})}{\varepsilon}\\
\end{split}\end{equation}


Решение уравнения сводиться к решению линейной системы:
\begin{equation*}\begin{split}
\bar\alpha&=\bar\alpha_{n+1}-\bar\alpha_{n}\\
X' \bar \alpha &=-\gamma X\\
\begin{pmatrix} a & b \\ c & d \end{pmatrix}\cdot\begin{pmatrix} \alpha_1 \\ \alpha_2 \end{pmatrix}&= -\gamma \begin{pmatrix} X_1 \\ X_2 \end{pmatrix}
\end{split}\end{equation*}
то есть:
\begin{equation*}\begin{split}
a\cdot \alpha_1 + b \cdot \alpha_2 = -\gamma \cdot X_1\\
c\cdot \alpha_1 + d \cdot \alpha_2 = -\gamma \cdot X_2
\end{split}\end{equation*}

а именно:
\begin{equation*}\begin{split}
\alpha_1 = -\gamma \frac{X_1 \cdot d - X_2 \cdot b}{a \cdot d - c \cdot b}\\
\alpha_2 = -\gamma \frac{X_1 \cdot c - X_2 \cdot a}{c \cdot b - a \cdot d}
\end{split}\end{equation*}


Следовательно следующий параметр пристрелки считается по формулам:
\begin{equation*}\begin{split}
\alpha_{{n+1}_1} &=\alpha_{{n}_1} -\gamma h_1\\
\alpha_{{n+1}_2} &=\alpha_{{n}_2} -\gamma h_2\\
h_1 &=\frac{X_1 \cdot d - X_2 \cdot b}{a \cdot d - c \cdot b} \\
h_2 &= \frac{X_1 \cdot c - X_2 \cdot a}{c \cdot b - a \cdot d}\\
\end{split}\end{equation*}

При выборе значения $\gamma$ сравниваем нормы функций невязок для $\bar \alpha_{n+1}$ и $\bar \alpha_n$.

В начале считаем $\gamma = 1$. Если не выполняется следующее неравенство: \begin{equation}\label{nev1}||X(\bar\alpha_{n+1})||<||X(\bar\alpha_n)||.\end{equation}
то $\gamma$ умножается на $\frac{1}{2}$ и параметры пристрелки пересчитываются для такого значения $\gamma$. В случае, если неравенстве
На первом шаге $\gamma=1$. Если после первой итерации неравенство \ref{nev1} вновь не выполнено, то $\gamma$ умножается на $\frac{1}{2}$ и далее, пока либо не подберем хорошие параметры пристрелки, либо $\gamma$ не станет меньше выбранной локальной погрешности $\varepsilon$.
\section{Тест решения задачи Коши - гармонический осциллятор.}
Проверим выбранный метод численного интегрирования при решении задачи с известным решением, а именно системы дифференциальныхИс уравнений гармонического осциллятора:  $\left\{ \begin{matrix} \dot x =& y,\\ \dot y =& -x. \end{matrix}\right.$  с начальными условиями  $\left\{ \begin{matrix}  x(0) =& 0,\\  y(0) =& 1. \end{matrix}\right.$

Используемый метод численного интегрирования - явный метод Рунге-Кутта с оценкой погрешности на шаге через 4 производную для различных значений $T$ и различных значений максимально допустимой относительной погрешности на шаге интегрирования $\varepsilon$. Интегрирование ведется на отрезке $[0,T]$. Функции невязок - $|x(T)|$ и $|y(T) - \cos{T}|$; число шагов интегрирования - $steps$; $\Delta x (t) , \Delta y (t)$ - максимальное отклонение полученных значений от аналитического решения  $\left\{ \begin{matrix} x =& \sin{t},\\ y =& \cos{t}. \end{matrix}\right.$, где максимум берется по всем шагам. $\delta_K(T)$ - оценка глобальной погрешности, которая считается по формуле $$\delta_{n+1}=E_n+\delta_n e^l, где l=\int\limits_{x_n}^{x_n+h}\lambda(v) \,dv$$
$$\delta_0=0$$
где $E_n$ - главный член в оценке локальной погрешности $\lambda(v)$ - логарифмическая норма матрицы Якоби $\mathcal{J}$ системы дифференциальных уравнений, равная наибольшему собственному числн матрицы $\frac{1}{2}(\mathcal{J}+\mathcal{J}^T)$.

Как было посчитано в пункте \ref{lognorm}, \ $\lambda(x) = 0$ и глобальная погрешность приближенно вычисляется по формуле
$$\delta_{n+1}=E_n+\delta_n $$. Числа Рунге $R_x= \left| \frac{x_{10^{-8}}(T)-x_{10^{-10}}(T)}{x_{10^{-10}}(T)-x_{10^{-12}}(T)}\right|$,  $R_y= \left| \frac{y_{10^{-8}}(T)-y_{10^{-10}}(T)}{y_{10^{-10}}(T)-y_{10^{-12}}(T)}\right|$. %которые должны быть примерно равны $100^{\frac{4}{5}}$


\begin{center}
\begin{tabular}[t]{|c||c|c||c|c||c|c|}
\hline
$\varepsilon= 10^{-8}$ & $|x(T)|$ & $|y(T)-1|$  & $\Delta x(\cdot)$ & $\Delta y(\cdot)$ & $\delta(T)$ & $steps$  \\
\hline
&0 & & & & & \\
\hline
&0 & & & & & \\
\hline
&0 & & & & & \\
\hline
&0 & & & & &\\
\hline
&0 & & & & & \\
\hline
\end{tabular}

\begin{tabular}[t]{|c||c|c||c|c||c|c|}
\hline
$\varepsilon= 10^{-10}$ & $|x(T)|$ & $|y(T)-1|$  & $\Delta x(\cdot)$ & $\Delta y(\cdot)$ & $\delta(T)$ & $steps$  \\
\hline
&0 & & & & & \\
\hline
&0 & & & & & \\
\hline
&0 & & & & & \\
\hline
&0 & & & & &\\
\hline
&0 & & & & & \\
\hline
\end{tabular}

\begin{tabular}[t]{|c||c|c||c|c||c|c|}
\hline
$\varepsilon= 10^{-12}$ & $|x(T)|$ & $|y(T)-1|$  & $\Delta x(\cdot)$ & $\Delta y(\cdot)$ & $\delta(T)$ & $steps$  \\
\hline
&0 & & & & & \\
\hline
&0 & & & & & \\
\hline
&0 & & & & & \\
\hline
&0 & & & & &\\
\hline
&0 & & & & & \\
\hline
\end{tabular}
\end{center}

 \begin{center} \begin{tabular}{|c||c|c|}
 \hline
$ T $& $R_x(T)$ & $R_y(T)$ \\
\hline
$10 \pi $&  8.03357 & 8.03281   \\
\hline
$10^2 \pi $& 8.01877 & 7.97836  \\
\hline
$10^3 \pi$& 8.02843 & 8.02410  \\
\hline
$10^4 \pi$  & 7.59795 & 7.50527  \\
 \hline
 \end{tabular}
 \end{center}

\section{Сравнение логарифмической нормы и максимального сингулярного числа.}\label{lognorm}
Матрица Якоби системы дифференциальных уравнений имеет вид:

\begin{equation*}
\mathcal J=\begin{pmatrix} 0 & 1 & 0 & 0 \\ 0 & -e^{-\alpha x}+\alpha x e^{-\alpha x}  & 0 & 1 \\ -2 \alpha p_y e^{-\alpha x}+\alpha^2 p_y x e^{-\alpha x}  & 0 & 0 & e^{-\alpha x}-\alpha x e^{-\alpha x}\\ 0& 0 & -1&0 \end{pmatrix} \end{equation*}

Для определения скорости распространения ошибки в оценках глобальной погрешности определяется логарифмическая норма матрицы $\lambda (\mathcal J)$ - максимальное собственное значение матрицы $(\mathcal J+\mathcal J^T)/2$ и норма матрицы $||\mathcal  J ||$ - максимальное сингулярное число матрицы $\mathcal J^T \mathcal J$.

\begin{equation*}
\mathcal J^T=\begin{pmatrix} 0 & 0 &  -2 \alpha p_y e^{-\alpha x}+\alpha^2 p_y x e^{-\alpha x} & 0 \\ 1 &  -e^{-\alpha x}+\alpha x e^{-\alpha x}& 0 &  0\\ 0 & 0 & 0 &-1\\ 0& 1 & e^{-\alpha x}-\alpha x e^{-\alpha x} &0 \end{pmatrix} \end{equation*}

Обозначим за $B = -2 \alpha p_y e^{-\alpha x}+\alpha^2 p_y x e^{-\alpha x}$, и за $A=-e^{-\alpha x}+\alpha x e^{-\alpha x}$. Тогда


\begin{equation*}
\frac{1}{2}(\mathcal J^T +\mathcal J)=\begin{pmatrix} 0 & \frac{1}{2} &  \frac{B}{2} & 0 \\ \frac{1}{2} &  A & 0& \frac{1}{2} \\ \frac{B}{2} & 0 & 0 &-\frac{A+1}{2}\\ 0& \frac{1}{2} & -\frac{A+1}{2}&0 \end{pmatrix} \end{equation*}

Поскольку нахождение логарифмической нормы матрицы в явном виде невозможно для данной системы, то воспользуемся следующим критерием: максимальное собственное число матрицы не превосходит максимального по модулю значения матрицы, то есть $$\lambda(A)\le \max{|a_{ij}|}$$

Тогда максимальное сингулярное число находим как максимальное собственное значение следующей матрицы:


\begin{equation*}
\mathcal J^T \mathcal J=\begin{pmatrix}  B^2 & 0 & 0 & -AB \\ 0 & 1+A^2 & 0 & A \\ 0 & 0 & 1 & 0\\ -AB& 0 &0 &1+A^2 \end{pmatrix}
\end{equation*}

\begin{equation*}
det(J^T J -x E)=(1-x)(-x^3+x^2(B^2+2+2A^2)+x(-2B^2-A^2B^2-1-A^4-A^2)+B^2)=0.
\end{equation*}

....................




\section{Аналитическое решение задачи.}\label{pissdeads}
Итак, вернемся к рассмотрению задачи \ref{system} и попробуем решить ее аналитически.
Дифференцируя первое уравнение системы, имеем

 $$\ddot x = \dot y = p_y - x e^{-\alpha }$$

$$\dddot x = \dot p_y - e^{-\alpha x} +\alpha x e^{-\alpha x} = -p_x - e^{-\alpha x} +\alpha x e^{-\alpha },$$

$$x^{(4)}=-\dot p_x +2 \alpha e^{-\alpha x}- \alpha^2 x e^{-\alpha x} = -p_y e^{-\alpha x} + \alpha p_y x e^{-\alpha x}+ 2 \alpha e^{-\alpha x}- \alpha^2 x e^{-\alpha x},$$

$$x^{(4)}+x^{(2)}(e^{-\alpha x} - \alpha x e^{-\alpha x}) - \alpha^2 x e^{-2\alpha x}+\alpha^2 x^2 e^{-\alpha x}+x e^{-2\alpha x}-2\alpha e^{-\alpha x} = 0.$$

Поскольку для произвольного значения $\alpha$ аналитически такое дифференциальное уравнение не решается, то решим его для значения $\alpha = 0$, и именно этот результат далее сравним с численным решением.

В этом случае система сводится к одному уравнению: $$x^{(4)}+2x^{(2)}+x=0$$

Характеристический полином: $$\lambda^4+2\lambda^2+1=(\lambda^2+1)^2=0.$$

Его корни имеют вид $$\lambda_{1,2,3,4}=\pm i$$

Тогда решение запишется в виде: $$x(t) = (A+Bt)\cos{t}+(C+Dt)\sin{t}$$

И решение системы имеет вид:
$$\left\{
\begin{aligned}
&x(t) =(A+Bt)\cos{t}+(C+Dt)\sin{t}\\
&y(t)= (B+C)\cos{t}+Dt\cos{t}+(D-A)\sin{t}-B t \sin{t},\\
&p_x(t) = 2 D\sin{t}+2 B \cos{t},\\
&p_y(t) = -2 B\sin{t}+2 D \cos{t}
\end{aligned}\right. $$

Из начальных и граничных условий имеем:

$$\left\{
\begin{aligned}
&x(\pi/2) =  C+D \frac{\pi}{2} = 0,\\
&y(0)= B + C= 0,\\
&y(\pi/2) = D - A - B\frac{\pi}{2}= 1,\\
&p_x(0) = x(0) = 2 B = A
\end{aligned}\right. $$


Тогда:

$$\left\{
\begin{aligned}&A = \frac{\pi}{1-\pi-\frac{\pi^2}{4}},\\
& B = \frac{\pi}{2(1-\pi-\frac{\pi^2}{4})},\\
&C = - \frac{\pi}{2(1-\pi-\frac{\pi^2}{4})},\\
&D = \frac{1}{1-\pi-\frac{\pi^2}{4}}
\end{aligned}\right. $$

     Итак, можем предоставить точное решение:

$$
\begin{cases}
&x(t) =(\frac{\pi}{1-\pi-\frac{\pi^2}{4}}+\frac{\pi}{2(1-\pi-\frac{\pi^2}{4})}t)\cos{t}+(- \frac{\pi}{2(1-\pi-\frac{\pi^2}{4})}+\frac{1}{1-\pi-\frac{\pi^2}{4}}t)\sin{t}\\
&y(t)= (\frac{\pi}{2(1-\pi-\frac{\pi^2}{4})}+- \frac{\pi}{2(1-\pi-\frac{\pi^2}{4})})\cos{t}+\frac{1}{1-\pi-\frac{\pi^2}{4}}t\cos{t}+(\frac{1}{1-\pi-\frac{\pi^2}{4}}-\frac{\pi}{1-\pi-\frac{\pi^2}{4}})\sin{t}-\frac{\pi}{2(1-\pi-\frac{\pi^2}{4})} t \sin{t},\\
&p_x(t) = 2 \frac{1}{1-\pi-\frac{\pi^2}{4}} \sin{t}+2 \frac{\pi}{2(1-\pi-\frac{\pi^2}{4})}\cos{t},\\
&p_y(t) = -2 \frac{\pi}{2(1-\pi-\frac{\pi^2}{4})} \sin{t}+2 \frac{1}{1-\pi-\frac{\pi^2}{4}} \cos{t}
\end{cases} $$

Точное значение параметров пристрелки: $x(0) = \frac{\pi}{1-\pi-\frac{\pi^2}{4}}\approx -0.681622241504944, p_y(0) =  \frac{\pi}{1-\pi-\frac{\pi^2}{4}}\approx -0.433934196227558$, точное значение интеграла $B_0 = \pi(( \frac{\pi}{2(1-\pi-\frac{\pi^2}{4})})^2+(\frac{1}{1-\pi-\frac{\pi^2}{4}})^2)-4 \frac{\pi}{2(1-\pi-\frac{\pi^2}{4})} \frac{1}{1-\pi-\frac{\pi^2}{4}}+ (\frac{\pi}{1-\pi-\frac{\pi^2}{4}})^2\approx 0.681622241504944$.
\section{Численное решение задачи.}\label{chis}
Задача была решена для каждого значения параметра $\alpha$. За максимально допустимую относительную погрешность на шаге решения $\varepsilon$ принималось $10^{-8}, 10^{-10}, 10^{-12}$.

В приведенной таблице даны результаты решения для различных значений $\alpha$. Здесь $x(0), p_y(0)$ - найденные методом стрельбы начальные условия, $|x(\pi/2)|, |y(\pi/2)-1|$ - функции невязки, $I$ - значение интеграла.
\begin{center}
\begin{tabular}[t]{|c|c||c|c||c|c|}
\hline
$\varepsilon= 10^{-6}$ & $\alpha$ & $x_0$ & $p_{y}(0)$ & $|x(\pi/2)|$ & $|y(\pi/2)-1|$  \\
\hline
&1.0 & -0.513631025928664 & -0.953204972729211& 6.20366e-07 &2.05637e-07  \\
\hline
&5.0& -0.055086724631734 &-1.224772565144717 &  4.63981e-09& 2.18205e-09\\
\hline
&10.0 & & & - &  \\
\hline
&15.0 & & & - &  \\
\hline
&20.0 & & & - &  \\
\hline
\end{tabular}

~\

~\
\begin{tabular}[t]{|c|c||c|c||c|c|}
\hline
$\varepsilon= 10^{-10}$ & $\alpha$ & $x_0$ & $p_{y}(0)$ & $|x(\pi/2)|$ & $|y(\pi/2)-1|$  \\
\hline
&1.0& -0.513630921789703 &-0.953205888206639& 4.85116e-09& 1.609e-09\\
\hline
&5.0 &-0.055086722032916&-1.224772561237645 & 7.23278e-11 &3.39089e-11  \\
\hline
&10.0 & & & - &  \\
\hline
&15.0 & & & - &  \\
\hline
&20.0 & & & - &  \\
\hline
\end{tabular}

~\

~\
\begin{tabular}[t]{|c|c||c|c||c|c|}
\hline
$\varepsilon= 10^{-12}$ & $\alpha$ & $x_0$ & $p_{y}(0)$ & $|x(\pi/2)|$ & $|y(\pi/2)-1|$  \\
\hline
&0.0& $1.2 \cdot 10^{-7}$ & & & \\
\hline
&5.0 &-0.055086721994115 & -1.224772561165411 &  5.70763e-13 &5.51115e-13  \\
\hline
&10.0 & -0.513630921154440 & -0.953205893840170 & 7.62393e-11& 2.53808e-11 \\
\hline
&15.0 & & & - &  \\
\hline
&20.0 & & & - &  \\
\hline
\end{tabular}

~\

~\
\begin{tabular}[t]{|c||c|c|c|}
\hline
 $\alpha$ & $I$ для $\varepsilon = 10^{-8}$ &  $I$ для $\varepsilon = 10^{-10}$ &  $I$ для $\varepsilon = 10^{-12}$ \\
\hline
1 & 0.723178547693589 & 0.723178538019265 & 0.723178537944255 \\
\hline
5. & 0.592484331346220  & 0.592484323394653 & 0.592484323333014\\
\hline
10.0 &  0.551142840794595  &0.551142833231005& 0.551142833172363 \\
\hline
15.0 &  0.535401571629420 &   0.535401564199650&  0.535401564142042 \\
\hline
22.0 &  0.524837925847553& 0.941341332291957& 0.524837925790687\\
\hline
\end{tabular}
\end{center}

В следующей таблице приведены значения глобальной погрешности в конце интегрирования для различных значений $\alpha$.

\begin{center}\begin{tabular}[t]{|c||c|c|c|}
\hline
$\alpha$ & $\delta(\pi/2)$ при $10^{-8}$ & $\delta(\pi/2)$ при $10^{-10}$& $\delta(\pi/2)$ при $10^{-12}$ \\
\hline
0 & & &  \\
\hline
5. & &  & \\
\hline
10.0 & & & \\
\hline
20.0  & & & \\
\hline
25.0 & & & \\
\hline
\end{tabular}
\end{center}
\section{Сравнение численного и аналитических решений.}
В пункте \ref{pissdeads} задача \ref{system} была решена аналитически при $\alpha = 0$. В приведенной таблице указаны отклонения параметров пристрелки и максимумы отклонения значений переменных $x, y, p_x, p_y$ для точек, в которые был сделан шаг интегрирования, при различных значениях максимально допустимую относительную погрешность на шаге решения $\varepsilon$ : $10^{-8}, 10^{-10}, 10^{-12}$.

\begin{center}

\begin{tabular}[t]{|c||c|c||c|c|c|c|}
\hline
 & $\Delta x(0)$ &$\Delta y(0)$ & $\Delta x(\cdot)$ & $\Delta y(\cdot)$ & $\Delta p_x(\cdot)$  & $\Delta p_y(\cdot)$  \\
\hline
$\varepsilon = 10^{-8}$ & & & & & & \\
\hline
$\varepsilon = 10^{-10}$ & & & & & & \\
\hline
$\varepsilon = 10^{-12}$  & & & & & & \\
\hline
\end{tabular}
\end{center}
\section{Оценка погрешности.}
Пусть $E_n$ - локальная погрешность на шаге $n$. Глобальная погрешность считается по следующей рекуррентной формуле $$\delta_{n+1}=E_n+\delta_n e^l, где l=\int\limits_{x_n}^{x_n+h}\lambda(v) \,dv$$
$$\delta_0=0$$
где $\lambda(v)$ это наибольшее собственное число матрицы $\frac{1}{2}(\mathcal{J}+\mathcal{J}^T)$. $\mathcal{J}$ - матрица Якоби системы дифференциальных уравнений.\

Поскольку в данной задаче получается лишь оценить $\lambda(v)$, то считаем $\lambda(v)\le \max \limits_{ij}|a_{ij}|$ , здесь $(a_{ij}) = \frac{1}{2}(\mathcal{J}+\mathcal{J}^T)$

Таким образом, $\delta_{n+1} \le E_n+\delta_n e^{h \max \limits_{ij}|a_{ij}|}$

\section{Правило Рунге.}
При различных значениях $\alpha$ выбирается 4 различных момента времени и вычисляется число Рунге для всех фазовых и сопряженых переменных (то есть для
$x,y, p_x, p_y$). Число Рунге считаем следующим образом:
$$R=\frac{X_{10^{-8}}-X_{10^{-10}}}{X_{10^{-10}}-X_{10^{-12}}}$$
 где $X_{10^{-8}}$ -  значение фазовой переменной $X$ при $\varepsilon = 10^{-8}$
 В таблице приведены числа Рунге для моментов времени $T/4, T/2, 3T/4, T$ (здесь $T$ - отрезок, на котором изменяется $t$, то есть $(0, \frac{\pi}{2})$) и оценка глобальной погрешности для значений $\alpha=0, 5, 10, 20, 25$.
 \begin{center} \begin{tabular}{|c|c||c|c|c|c|}
 \hline
$\alpha$ & $X$ & $R_{T/4}$ & $R_{T/2}$ & $R_{3T/4}$ & $R_T$ \\
\hline
0.1 & $x$ & 8.03357 & 8.03281 & 15.65843 & 8.24146  \\
\hline
0.1 & $ p_x $ & 8.01877 & 7.97836 & 18.12140& 8.06067  \\
\hline
0.5 & $y$ & 8.02843 & 8.02410 & 15.64284 & 8.58319 \\
\hline
0.5 & $p_y $ & 7.59795 & 7.50527 & 58.81771 & 8.57883  \\
 \hline
 \end{tabular}
 \end{center}


\section{Просчет назад.}
В следующих таблицах приведены результаты "просчёта назад". Решалась задача Коши на промежутке от $0$ до $\pi/2$. В качестве начальных данных берем числа, указанные в условии задачи, и числа, полученные в результате пристрелки. Далее полученные значения переменных в точке $\pi/2$  принимались в качестве начальных условий для интегрирования аналогичных уравнений от $\pi/2$ до $0$.

Далее $x(0), p_y(0)$ - значения, полученные в результате метода стрельбы (при $\varepsilon = 10^{-12}$), $y(0), p_x(0)$ взяты из условия задачи, $x'(0), y'(0), p_x'(0), p_y'(0)$ - значения переменных в нуле, полученные в результате просчета назад. Максимально допустимая относительная погрешность на шаге решения - $\varepsilon = 10^{-12}$

\begin{center}
\begin{tabular}[t]{|c||c|c|}
\hline
$\alpha$ & $x(0)$ & $x'(0)$  \\
\hline
1.0 & & \\
\hline
5.0 &  &   \\
\hline
10.0 &  &  \\
\hline
15.0 & &   \\
\hline
22.0 & &\\
\hline
\end{tabular}

~\

~\
\begin{tabular}[t]{|c||c|c|}
\hline
$\alpha$ & $y(0)$ & $y'(0)$  \\
\hline
1.0 & 0.0 & \\
\hline
5.0 & 0.0 &  \\
\hline
10.0 & 0.0 &    \\
\hline
15.0 & 0.0 & \\
\hline
22.0 & 0.0 &   \\
\hline
\end{tabular}

~\

~\
\begin{tabular}[t]{|c||c|c|}
\hline
$\alpha$ & $p_x(0)$ & $p_x'(0)$  \\
\hline
1.0 &  & \\
\hline
5.0 &  &  \\
\hline
10.0 & &  \\
\hline
15.0 & &   \\
\hline
22.0 & &  \\
\hline
\end{tabular}

~\

~\
\begin{tabular}[t]{|c||c|c|}
\hline
$\alpha$ & $p_y(0)$ & $p_y'(0)$  \\
\hline
1.0 & & \\
\hline
5.0 & & \\
\hline
10.0 & &  \\
\hline
15.0 & &   \\
\hline
22.0 & & \\
\hline
\end{tabular}
\end{center}

Нетрудно заметить, что различия имеют порядок примерно $10^{-13}$.

\section{Графики.}
ГРАФИКИ
\section{Исследование оптимальности экстремалей.}

Исследуем оптимальности в задаче Лагранжа \ref{lag} в соответствии c алгоритмом, подробно рассмотренном в [1]. Рассмотрим пространство управляемых процессов $X = C^1([0,\frac{1}{2}], \mathbf R) \times C([0, \frac{\pi}{2}])$, в котором рассматривается задача. Его элементы называем управляемыми процессами и обозначаем $\xi :=\{x(\cdot), y(\cdot), u(\cdot)\}$.
Обозначим $\hat \xi:=\{\hat x(\cdot),\hat y(\cdot),\hat u(\cdot)\}$ элемент этого пространства, который удовлетворяет необходимым и достаточным условиям оптимальности второго порядка и находится в результате решения краевой задачи.

Оператор, действующий из пространства $X$ в пространство $Y = C([0, \frac{\pi}{2}], \mathbf R^2) \times \mathbf R^2$, имеет следующий вид

\begin{equation}\label{system2}
\left\{
\begin{aligned}
&\dot x = y,\\
&\dot y = u - x \exp(-\alpha x),\\
&y(0) = 0,\\
&y(\pi/2) = 1,\\
&x(\pi/2)=0,\\
&u \in \mathbf R,\\
&\alpha \in [0.0, 25.0]
\end{aligned}\right.
\end{equation}

Рассмотрим многообразие, определяемое в точке $\hat \xi$ пространства $X$ условиями  \ref{system2}, построим к нему касательное пространство и его элементы обозначим $\delta \xi = \{\delta x, \delta y, \delta u \}$.

В силу теоремы Люстернака такое пространство совпадает с ядром оператора \ref{system2}. Это накладывает на $\delta \xi$ следующие условие
\begin{equation}\label{sys}
\left\{
\begin{aligned}
&\dot \delta x = \delta y,\\
&\dot \delta y = \delta u +(-e^{-\alpha x}+\alpha x e^{-\alpha x})\delta x,\\
&\delta y(0) = 0,\\
&\delta y(\pi/2) = 1,\\
&\delta x(\pi/2)=0,\\
&\delta u \in \mathbf R,\\
&\alpha \in [0.0, 25.0]
\end{aligned}\right.
\end{equation}

Вторая вариация функции Лагранжа имеет вид \begin{equation}\label{L}\mathcal L_{\xi \xi} (\delta \xi, \delta \xi) = \int \limits_0^{\frac{\pi}{2}}(\hat p_y(-2 \alpha e^{-\alpha \hat x}+ \alpha^2 \hat x e^{-\alpha \hat x})(\delta x)^2 + (\delta u)^2)dt\end{equation}

Тогда необходимое условие оптимальности перепишется в виде $\mathcal L_{\xi \xi}[\delta \xi, \delta \xi]\ge 0$, а достаточное - в виде $\mathcal L_{\xi \xi}[\delta \xi, \delta \xi]\ge \varepsilon \varphi(\delta \xi)$, где $\varphi(\delta \xi) = \int \limits_0^{\frac{\pi}{2}} (\delta u)^2 dt$. Отметим, в общем случае $\varphi(\delta \xi)$ содержит еще и член вида $x^2(\tau)$, где $\tau$ - любая точка из $[0, \frac{\pi}{2}]$. Поскольку временной отрезок фиксирован, то считаем этот член нулевым.

Для неотрицательности $\mathcal L_{\xi \xi}[\delta \xi, \delta \xi]$ необходимым является условие Лежандра:
$\mathcal L_{uu}[v,v]\ge 0 \qquad \forall v$. В данной задаче оно выполняется в усиленной форме: $\mathcal L_{uu}[v,v] = 1$.
Поэтому для неотрицательности второй вариации необходимо отсутствие на интервале $(0,\frac{\pi}{2})$ сопряженных точек (условие Якоби),
для положительности второй вариации достаточно отсутствие на полуинтервале $(0,\frac{\pi}{2}]$ сопряженных точек (усиленное условие Якоби)

Сопряженной точкой называем $\tau$ - верхнюю грань таких $t_{*}>t_0$, что вторая вариация Лагранжа на любых двух элементах касательного пространства строго положительна, а на отрезке $[t_{*}, t_1]$ все $\delta \xi \equiv 0$.

Такая точка определяется следующей системой уравнений:

\begin{equation}\
\left\{
\begin{aligned}
&\dot \delta x = \delta y,\\
&\dot \delta y = \delta u +(-e^{-\alpha \hat x}+\alpha x e^{-\alpha \hat x})\delta x,\\
&\dot q_x = q_y(e^{-\alpha \hat x}+\alpha \hat x e^{-\alpha \hat x})-(\delta x)p_y(2\alpha e^{-\alpha x}-\alpha^2 x e^{-\alpha x},\\
&\dot q_y = - q_x,\\
&\delta u = q_y,\\
&\delta y(\tau) = 1,\\
&\delta y(0) = 0,\\
&\delta x(\tau)=0,\\
&q_x(0) = 0,\\
&\delta u \in \mathbf R,\\
&\alpha \in [0.0, 25.0]
\end{aligned}\right.
\end{equation}

Эта система может быть формально получена в результате применения необходимых условий оптимизации к задаче минимизации функционала (\ref{L}) при выполнении системы линейных условий (\ref{sys}) на временном интервале $[0,\tau]$. Множитель Лагранжа был выбран $\frac{1}{2}$.

Так как система дифференциальных уравнений рассматриваемой краевой задачи линейна и однородна, то все решения задачи Коши, удовлетворяющие условиям при $t=0$ ($\delta y (0)=0, q_x(0)=0)$ могут быть получены  в виде линейной комбинации двух линейно независимых решений задач Коши с начальными условиями:
\begin{equation}\label{2}\begin{pmatrix} \delta x^1(0)\\ \delta y^1(0)\\ q^1_x(0) \\ q^1_y(0) \end{pmatrix}=\begin{pmatrix} 1\\ 0 \\ 0 \\ 0 \end{pmatrix} \qquad \begin{pmatrix} \delta x^2(0)\\ \delta y^2(0)\\ q^2_x(0) \\ q^2_y(0) \end{pmatrix}=\begin{pmatrix} 0\\ 0 \\ 0 \\ 1 \end{pmatrix} \end{equation}.

Если в некоторый момент времени $\tau \in (0, \frac{\pi}{2})$ выполняется следующее соотношение
\begin{equation}\label{det}det
\begin{pmatrix}\delta x^1(\tau) & \delta x^2(\tau)\\\delta y^1(\tau) & \delta y^2(\tau)\end{pmatrix} = 0
\end{equation}
то краевая задача имеет нетривиальное решение, а точка $\tau$ будет сопряженной.

Поскольку матрица Якоби данной системы совпадает с матрицей Якоби краевой задачи (\ref{system}), ее решение совпадает с решением краевой задачи.

Как получили в пункте \ref{pissdeads}, система имеет следующее аналитическое решение для $\alpha = 0$:

$$\left\{
\begin{aligned}
&\delta x(t) = (A'+B't)\cos{t}+(C'+D't)\sin{t},\\
&\delta y(t)= (B+C)\cos{t}+Dt\cos{t}+(D-A)\sin{t}-B t \sin{t},\\
&\delta p_x(t) = 2 D\sin{t}+2 B \cos{t},\\
&\delta p_y(t) = -2 B\sin{t}+2 D \cos{t}
\end{aligned}\right. $$


Таким образом, начальные условия (\ref{2}) дают следующие соотношения для констант в решении:

$$\left\{
\begin{aligned}
&A'_1 = 1,\\
&B'_1+C'_1 = 0,\\
& 2 B'_1=0,\\
&2 D'_1=0
\end{aligned}\right. $$

$$\left\{
\begin{aligned}
&A'_2 = 0,\\
&B'_2+C'_2 = 0,\\
&2 B_2'=1,\\
&2 D_2' =0
\end{aligned}\right. $$

Можно найти константы:
$$\left\{
\begin{aligned}
&A'_1 = 1,\\
&B_1' = 0,\\
&C_1' = 0,\\
&D_1' = 0
\end{aligned}\right. $$

$$\left\{
\begin{aligned}
&A'_2 = 0,\\
&B_2' = \frac{1}{2},\\
&C_2' = -\frac{1}{2},\\
&D_2' = 0
\end{aligned}\right. $$

Тогда уравнение \ref{det} перепишется в виде $$- \tau\cos{\tau}\sin{\tau} \frac{1}{2}+  \tau \cos{\tau}\sin{\tau} \frac{1}{2}-\frac{1}{2}\sin{\tau}= 0$$


Решение этого уравнения: $\tau* = 0, \pi$.
Таким образом, на полуинтервале $(0, \pi/2]$ сопряженные точки отсутствуют, следовательно, на экстремали задачи выполняется усиленное условия Якоби и достигается слабый локальный минимум. Поскольку выполняется условие квазирегулярности (усиленное условие Вейерштрасса , функция $L$ является выпуклой функцией управления $u$, то есть $L_uu$ положительно определена в окрестности экстремали), то этот минимум является сильным локальным минимумом. Так как задача линейно-квадратичная, то сильный локальный минимум совпадает с абсолютным. Таким образом, для $\alpha = 0$ полученные экстремали оптимальны и доставляют минимум функции $B_0$.

 \newpage
{\huge \bf Литература.}
 \begin{itemize}
\item[1.] И. С. Григорьев. Методическое пособие по численным методам решения краевых задач принципа максимума в задачах оптимального управления.
\item[2.] В. В. Александров, Н. С. Бахвалов, К. Г. Григорьев, Г. Ю. Данков, М. И. Зеликин, С. Я. Ищенко, С. В. Конягин, Е. А. Лапшин, Д. А. Силаев, В. М. Тихомиров, А. В. Фурсиков. Практикум по численным методам в задачах оптимального управления.


    \item[3.] Э. Хайрер, С. Нёрсетт, Г. Ваннер. Решение обыкновенных дифференциальных уравнений.
     \item[4.] Исаев В. К., Сонин В. В. Об одной модификации метода Ньютона численного решения краевых задач.

\end{itemize}
%\section{Оценка... 19 пункт}
\end{document}