\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[utf8]{inputenc}

\usepackage[russian]{babel}
\usepackage[intlimits]{amsmath}
\interdisplaylinepenalty=2500
\usepackage{amsfonts}
\usepackage{amsthm}
\usepackage{amssymb}
\usepackage{euscript}
\usepackage{amsmath}
\usepackage{latexsym}
\usepackage{graphicx}
\newtheorem{theorem}{Теорема}
\newtheorem{proposition}{Предложение}
\newtheorem{statement}{Утверждение}
\newtheorem{corollary}{Следствие}
\renewcommand{\proofname}{Доказательство}
\newtheorem{definition}{Определение}

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

\begin{document}
\begin{titlepage}
    \begin{center}

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

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

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



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

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


            Вавилкин Андрей
        \vfill
        4 курс, группа 422
    \end{center}
    \vfill



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

\end{center}



%----------------------------


\section{Постановка задачи.}
Найдем зависимость $\beta$ от $\alpha$, где $\beta, \alpha$ удовлетворяют задаче Коши следующего вида
$$
\begin{cases}
\begin{cases}
\dfrac {dx}{dt}=\alpha \left( x^{3}-1\right) -y, \\
\dfrac {dy}{dt}=x-\alpha \left( y^{3}+1\right), \\
\end{cases}\\
0<t<2\beta,\\
t=0: x=0,y=-1,\\
t=\beta: x+y=0;\\
0<\alpha<1.
\end{cases}
$$

\section{Обзор методов решения.}
\subsection{Метод Рунге-Кутта.}
Интегрирование проводилось при помощи метода Рунге-Кутта 4 порядка (подробное описание можно найти в $[1]$).

Дадим краткий обзор данного метод для общего случая.

Пусть имеется задача Коши
$$ \left\{
\begin{array}{rcl}
y'=f(x,y(x)),\\
y(0)=y_0.\\

\end{array}
\right.$$
$f$ считаем достаточно гладкой функцией. Приближенное решение $y_{n+1}$ в узле $x_n+h$ представляется в виде
$$y_{n+1}=y_{n}+\frac{1}{6}(k_1+2k_2+2k_3+k_4),$$
где коэффицинты $k_1, k_2, k_3, k_4$ считаем по формулам
$$k_1=h f(x_n, y_n),$$
$$ k_2 = h f(x_n +\frac{1}{2}h, y_n + \frac{1}{2}k_1),$$
$$ k_3 = h f(x_n +\frac{1}{2}h, y_n + \frac{1}{2}k_2),$$
$$ k_4 = h f(x_n +h, y_n + k_3).$$

Контрольный член погрешности, известный как контрольный член Егорова, определяется по формуле $$E=\frac{2}{3}(k_1-k_2-k_3+k_4)$$
и имеет порядок $O(h^3)$

В данном методе фигурирует параметр $h$ - шаг интегрирования. В следующем параграфе опишем метод выбора подходящего значения $h$.
\subsection{Процедура автоматического выбора шага.}
%Пусть на $n$ шаге алгоритма имеем решение $y_n$ в точке $x_n$ и шаг $h$. Считаем $y_{n+1}$ в точке $x_{n+1}=x_n+h$. Далее, по формуле $(1)$ строим число $E$ и сравниваем его с наперед заданным $varepsilon$ (считаем $\varepsilon = 10^{-5}, 10^{-7}, 10^{-9}$). Имеем 3 случая
%\begin{itemize}
%\item[--] если $|E|>\varepsilon$, то отменяем $n+1$ шаг и делим $h$ на $2$, то есть далее считаем $y_{n+1}$ в точке $x_{n+1}=x_n+\frac{h}{2}$
%\item[--] если $|E|<\frac{\varepsilon}{32}$, то отменяем $n+1$ шаг и умножаем $h$ на $2$, то есть далее считаем $y_{n+1}$ в точке $x_{n+1}=x_n+2р$
%\item[--] если $\frac{\varepsilon}{32}<|E|<\varepsilon$, то $h$ менять не нужно и можем продолжать интегрирование.
%\end{itemize}

 Пусть на $n$ шаге алгоритма имеем решение $y_n$ в точке $x_n$ и шаг $h$. Считаем $y_{n+1}$ в точке $x_{n+1}=x_n+h$ двумя спосоами - первый раз делаем один шаг длины $h$ и получаем значение $y_{n+1}^1$, второй раз - два шага длины $\frac{h}{2}$ и получаем $y_{n+1}^2$. Считаем модуль разности этих двух значений $E=|y_{n+1}^1-y_{n+1}^2|$. Если он оказыввается меньше наперед заданного числа $\varepsilon$ (считаем $\varepsilon = 10^{-5}, 10^{-7}, 10^{-9})$, то отменяем шаг, берем новое значение $h$, равное старому значению $h$ деленному на 2, и продолжаем интегрирование до тех пор, пока $E$ не примет допустимое значение.

Аналогично, считаем $y_{n+1}$ $x_{n+1}=x_n+2h$ двумя способами и получаем 2 значения $y_{n+1}^1$ и $y_{n+1}^2$, первое из которых отвечаеь одному шагу длины $2h$, а второе - двум шагам длины $h$ каждый. Тогда если $E=|y_{n+1}^1-y_{n+1}^2|$ оказывается меньше $\frac{\varepsilon}{2}$, то шаг увеличиваем в 2 раза и продолжаем интегрирование, начиная с точки $x_n$, до тех пор, пока $E$ не пример допустимое значение.

\subsection{Оценка глобальной погрешности.}
Пусть $E_n$ - локальная погрешность на шаге $n$. Глобальная погрешность счиатется по следующей формуле $$\delta_{n+1}=E_n+\delta_n e^l,\qquad \mbox{где } l=\int\limits_{t_n}^{t_n+h}\lambda(v) \,dv$$
$$\delta_0=0$$
где $\lambda$ это наибольшее собственное число матрицы $ \frac{1}{2}(\mathcal{J}+\mathcal{J}^T) $. Где $\mathcal{J}$ - матрица Якоби системы дифференциальных уравнений. Нетрудно посчитать, что $\lambda= 3\alpha x^2$
 Тогда в данной задаче $l = 3h\alpha x^2 $. Глобальную погрешность можем приближенно вычислять по следующией формуле   $$\delta_{n+1}=E_n+\delta_n e^{3h\alpha x^2}$$


\section{Результаты.}
Задача была запрогроммирована на языке $C$ при помощи приложения $CodeBlocks$ для $Windows$, и была получена таблица значений параметра  $\beta$ для различных $\alpha$, лежащих в заданном диапазоне $(0,1)$.  По такой таблице - набору точек - были построен график данной зависимости в программе Gnuplot.
 %КАРТИНКИ

\section{Тестирование.}
При фиксированном $\alpha$ выбирается 4 различных момента времени и вычисляется число Рунге для переменной $x$). Число Рунге считаем следующим образом:
$$R=\frac{x_{10^{-4}}-x_{10^{-6}}}{x_{10^{-6}}-x_{10^{-8}}}$$
 где $x_{10^{-8}}$ -  значение $x$ при $\varepsilon = 10^{-8}$.\\
 В таблице приведены числа Рунге  и оценка глобальной погрешности для значений времени, когда $x+y=0$.
 \begin{center}
 \begin{tabular}{|c|c|c|c|c|c|}
 \hline
$\alpha$ & $t_{x+y=0}$ & $R_{x}$ & $R_{y}$ & $\delta_x$ & $\delta_y$  \\
\hline
0.03 & 0.815  & 9.5609 & 8.57358 & 0.001204397430514 & 0.000485337989932 \\
\hline
 0.27 & 1.075  & 11.20773 & 9.06763 & 0.001279028082222 & 0.000359953963765 \\
 \hline
  0.75 & 4.365  & 7.26068 & 6.95665 & 0.000792562628044 & 0.000450364111816 \\
 \hline
 \end{tabular}
 \end{center}

 Ниже на графике изображена искомая зависимость $\beta$ от $\alpha$.

 \begin{center}
\includegraphics[scale=0.3]{img}
 \end{center}

 \section{Проверка работы программы.}
 Рассмотрим задачу следующего вида
$$
\begin{cases}
\frac{d x}{dt}=y,\\
\frac{d y}{dt}=-x,\\
x(0)=0,\\
y(0)=1.
\end{cases}
$$
Она интегрируется, и ее решение $x = \sin{t}, y=\cos{t}$. Графиком решения являются замкнутые окружности.
Интегрируя ее вышеописанным методом при $0<t<100 \pi $ и 3 значениях $\varepsilon = 10^{-4}, 10^{-6}, 10^{-8}$ получим график, представленный на рисунке. При увеличении точности (то есть при уменьшеии $\varepsilon$) график становит более точным и представляет собой спираль с достаточно близко расположенными дугами.


 \newpage
{\huge \bf Литература.}
 \begin{itemize}
\item[1.] О. Б. Арушанянб, С. Ф. Залёткин. Решение систем обыкновенных дифференциальных уравнений методами Рунге-Кутта, 2014
\item[2.] В. Ф. Дьяченко. Вычислительный практикум


    \item[3.] Э. Хайрер, С. Нёрсетт, Г. Ваннер. Решение обыкновенных дифференциальных уравнений.
     \item[4.] Н. С. Бахвалов, Н.П. Жидков, Г. М. Кобельков. Численные методы

\end{itemize}

\end{document}