doc/report.tex

\documentclass{article}
\usepackage[utf8x]{inputenc}
\usepackage[english,russian]{babel}
\usepackage{mflogo}
\usepackage{listings}

\lstdefinelanguage[guile]{Scheme}
        {morekeywords={define,let,map,lambda,if,cond,or,and,else},
         morecomment=[l]{;},
         morestring=[b]{"}}
\lstset{mathescape=true}
\lstset{language=[guile]Scheme}
%\lstset{numbers=left,stepnumber=10,numberblanklines=false}
\lstset{basicstyle=\footnotesize}
\lstset{frame=tb,frameround=ffff}
\lstset{breaklines=true,breakautoindent=true}

\newcommand{\filename}[1]{\texttt{#1}}
\newcommand{\procname}[1]{\texttt{#1}}
\newcommand{\includeplot}[2]{\begin{figure}[hb]
    \centering
    \includegraphics{#1__#2-plot.mps}
    \caption{График $u(x)$ для функции преломления \eqref{#2-initial-data}}
\end{figure}}
\renewcommand{\epsilon}{\varepsilon}
\renewcommand{\phi}{\varphi}
\providecommand{\abs}[1]{\left \lvert{#1}\right \rvert}
\providecommand{\op}[1]{\mathrm{#1}}
\providecommand{\comp}{\circ}

\usepackage{amsmath,amssymb}

\usepackage{graphics}

\DeclareGraphicsExtensions{.mps,.pdf,.eps}
\DeclareGraphicsRule{.scm-plot.mps}{mps}{*}{}
\DeclareGraphicsRule{.scm-all-plots.mps}{mps}{*}{}

\numberwithin{equation}{section}

\begin{document}

\author{Дмитрий Джус}
\title{Курсовая работа по теме\\ «Обыкновенные дифференциальные уравнения»}
\maketitle
\thispagestyle{empty}

\clearpage
\tableofcontents

\clearpage
\part{Постановка задачи}

\section{Математический аспект}

Настоящая курсовая работа посвящена построению приближённого решения
краевой задачи для дифференциального уравнения вида

\begin{equation}\label{deq}
  \frac{d^2 u}{dx^2} + n(x) u(x) = 0
\end{equation}

По условию, $n(x)$ представляет собой кусочно-непрерывную функцию,
равную константе $k^2$ на интервалах $(-\infty; 0)$ и $(a; +\infty)$ и
непрерывной известной функции на $[0; a]$. На границах $[0; a]$
решение и его первая производная должны удовлетворять условиям
сшивания:

\begin{subequations}\label{conds}
  \begin{align}
    u(0-0)& = u(0+0)&
    u'(0-0)& = u'(0+0)&\label{conds-left} \\
    u(a-0)& = u(a+0)&
    u'(a-0)& = u'(a+0)&\label{conds-right}
  \end{align}
\end{subequations}

\section{Физический аспект}

Задача соответствует возбуждения плоского слоя неоднородной
немагнитной среды электромагнитным полем, вектор электрического поля
распространяется вдоль оси $x$. $k$ имеет смысл волнового числа в
однородной части пространства, $n(x)$ — показатель преломления
неоднородной среды.

Из физических соображений, решение справа от слоя $[0; a]$ должно
иметь вид $u(x) = Be^{ikx}$, а слева — $u(x) = e^{ikx} +
Ae^{-ikx}$. Задача состоит в приближённом построении решения внутри
неоднородного слоя и отыскании комплексных констант $A$ и $B$
(называемых, соответственно, коэффициентами отражения и прохождения
электромагнитного поля).

Действительная область значений функции $n(x)$ соответствует
отсутствию поглощения энергии в среде, так что критерием правильности
полученных приближённых значений служит энергетическое тождество:

\begin{equation}\label{energy}
  \abs{A}^2 + \abs{B}^2 = 1
\end{equation}

\clearpage
\part{Решение}

\section{Общие замечания}

Предлагаемые к реализации в курсовой работае приближённые методы были
описаны на алгоритмическом языке Scheme (простом и строго
функциональном диалекте Лиспа).

\section{Решение методом построения фундаментальной матрицы}

\subsection{Описание}

Дифференциальное уравнение второго порядка \eqref{deq} можно представить
в матричной форме, введя $v(x) = \frac{du}{dx}$:

\begin{equation}\label{mdeq}
  \frac{d}{dx}
  \begin{pmatrix}
    u(x) \\
    v(x)
  \end{pmatrix}
  = \begin{pmatrix}
    0 & 1 \\
    -n(x) & 0
  \end{pmatrix}
  \begin{pmatrix}
    u(x) \\
    v(x)
  \end{pmatrix}
\end{equation}

Решение системы \eqref{mdeq} с матрицей $A(x) = \bigl( \begin{smallmatrix} 0
  & 1 \\ -n(x) & 0 \end{smallmatrix} \bigr)$ можно найти в общем виде, если
построена фундаментальная матрица $\Omega(x)$:

\begin{equation}\label{fm}
\Omega(x) = \begin{pmatrix}
  \omega_{11}(x) & \omega_{12}(x) \\
  \omega_{21}(x) & \omega_{22}(x)
\end{pmatrix}
\end{equation}

$\Omega(x)$ является решением следующей матричной задачи Коши:

\begin{equation}\label{fmeq}
  \frac{d}{dx} \Omega(x) = A(x) \Omega(x)
\end{equation}

с условием $\Omega(0) = \bigl( \begin{smallmatrix} 1& 0\\ 0&
  1\end{smallmatrix} \bigr)$.

Решение \eqref{mdeq} для $\forall x \in [0; a]$ записывается в виде

\begin{equation}\label{fmsolution}
\begin{pmatrix}
  u(x) \\
  v(x)
\end{pmatrix} = 
\begin{pmatrix}
  \omega_{11}(x) & \omega_{12}(x) \\
  \omega_{21}(x) & \omega_{22}(x)
\end{pmatrix}
\begin{pmatrix}
  1 + A \\
  ik(1 - A)
\end{pmatrix}
\end{equation}

Использование краевых условий \eqref{conds} позволяет получить
следующую линейную алгебраическую систему относительно неизвестных $A$
и $B$:

\begin{equation}\label{abeq}
  \begin{split}
    (\omega_{11}(a) - \omega_{12}(a) ik) A - e^{ika} B &= -ik
    \omega_{12}(a) - \omega_{11}(a) \\
    (\omega_{21}(a) - \omega_{22}(a) ik) A - ike^{ika} B &= -ik
    \omega_{22}(a) - \omega_{21}(a)
  \end{split}
\end{equation}

Нахождение значений компонентов матрицы $\Omega$ на правом конце $[0;
a]$, таким образом, позволит получить искомые значения $A$ и $B$.

Для поиска фундаментальной матрицы в точках $[0; a]$ предлагается
использование следующего численного метода, основанного на приближении
функции $n(x)$ кусочно-постоянной, то есть замене неоднородной среды
на плоско-слоистую с большим числом однородных слоёв.

Рассматриваемый интервал разбивается на $N$ интервалов точками
$x_i = \frac{a}{N}i$, $i = 1 \dotso N$, в середине каждого из
которых берётся точка $\overline{x_i} = \frac{a}{N}(i+\frac{i}{2})$. В
данной точке функция $n(x)$ аппроксимируется ступенчатой функцией
$\overline{n}(x) = n(\overline{x_i})$, так что \eqref{mdeq} обращается
в систему с постоянной матрицей $A_i = \bigl(
\begin{smallmatrix}0& 1 \\ -\overline{n}(x)& 0 \end{smallmatrix}
\bigr)$. Тогда на каждом интервале $[x_i; x_{i+1}]$ решение
\eqref{fmeq} может быть найдено явно в виде матричной экпоненты:

\[
\Omega_i(x) = e^{A_i (x - x_i)}
\]

При этом, значения фундаментальной матрицы $\Omega_N(x_{i+1})$ на
правом конце $[x_i; x_{i+1}]$ служат начальным условием задачи Коши на
$[x_{i+1}; x_{i+2}]$, так что для каждой точки $\hat{x} \in [x_i;
x_{i+1}]$ фундаментальная матрица \eqref{fm} имеет вид:

\begin{equation}\label{fmx}
  \Omega_N(\hat{x}) = e^{A_i(\hat{x}-x_i)} e^{A_{i-1}(x_i-x_{i-1})} \dotsm e^{A_0(x_1-x_0)}
\end{equation}

На правом конце интервала $[0; a]$ приближённая фундаментальная
матрица записывается в виде:

\begin{equation}\label{fma}
  \Omega_N(a) = e^{A_{N-1}(x_N-x_{N-1})} e^{A_{N-2}(x_{N-1}-x_{N-2})} \dotsm e^{A_0(x_1-x_0)}
\end{equation}

Для нахождения матричной экспоненты используется формула Тейлора:

\begin{equation}\label{matrix-exp}
\begin{split}
  e^{A_i(x-x_i)} \approx E +& \frac{1}{1!}{A_i(x-x_i)} +
  \frac{1}{2!}{A_i^2(x-x_i)^2}\\
  +& \frac{1}{3!}{A_i^3(x-x_i)^3} + \frac{1}{4!}{A_i^4(x-x_i)^4}
\end{split}
\end{equation}

После построения фундаментальной матрицы в $a$ коэффициенты $A$ и $B$
находятся из \eqref{abeq} с использованием полученных компонент
$\omega_{ij}^N(a)$.

Решение на каждом отрезке разбиения $[0; a]$ находится из
\eqref{fmsolution} с использование найденного $A$ и последовательности
приближений фундаментальной матрицы на $[x_i; x_{i+1}]$ для $\forall
i$.

\clearpage
\subsection{Реализация}

Описание основных процедур, необходимых для реализации метода,
находится в файле \filename{fundmatrix-solution.scm} (его полное содержание
приведено на странице \pageref{fundmatrix-solution.scm-full-listing}).

В исходном тексте программы активно используются последовательные
определения функций в терминах друг друга.

Так, основной вычислительный процесс метода — построение приближений
фундаментальной матрицы — выражен в функции
\procname{build-fundamentals}:

\input{build-fundamentals__fundmatrix-solution.scm-tag-listing}

\procname{build-fundamentals} представляет последовательность вычислений
фундаментальной матрицы на интервалах $[x_i; x_{i+1}]$ в виде частного
случая функции \procname{evolve-sequence}, возвращающей
последовательность $a_{s_1}, \dotsc, a_{s_n}$, где $a_{s_k} =
f(a_{s_{k-1}}, s_{k-1})$
(см. с. \pageref{shared.scm-full-listing}). Из определения
\procname{build-fundamentals} видно, что роль $s_k$ выполняют точки $x_i
= \frac{a}{N}i$ с $i = 1 \dotso N$, а $a_{s_{k-1}}$ — приближение
фундаментальной матрицы на предыдущем отрезке (на первом шаге
$\Omega_N$ приближается единичной матрицей). Комбинируются же эти
значения при помощи матричного умножения, как следует из \eqref{fmx}.

Функция \procname{variable-matrix} возвращает матрицу $A$ системы
\eqref{mdeq} для заданной по условию функции $n(x)$:

\input{variable-matrix__fundmatrix-solution.scm-tag-listing}

Возвращённая матрица потом используется в \procname{build-fundamentals}.

В \procname{build-fundamentals} используется и функция
\procname{matrix-exp}, вычисляющая матричную
экспоненту $f(x, y, z) = e^{A(x)(y-z)}$ путём разложения в ряд Тейлора
\eqref{matrix-exp} с использованием схемы Горнера:

\input{matrix-exp__matrices.scm-tag-listing}

В этой функции последовательность коэффициентов матричного многочлена
$\frac{1}{0!}, \frac{1}{1!}, \frac{1}{2!}, \dotsc$ генерируется
функцией \procname{exp-series-coefficients},
выраженной в терминах вышеупомянутой \procname{evolve-sequence}.

Вычисление матричной экспоненты в каждой точке выражено через
функцию высшего порядка \procname{general-horner-eval}, представляющую
собой обобщённую схему Горнера:

\input{general-horner-eval__shared.scm-tag-listing}

\procname{general-horner-eval}, в свою очередь, является частным случаем
функции высшего порядка \procname{fold-right}, реализующей свёртку
последовательности справа налево.

После построения последовательности фундаментальных матриц
осуществляется поиск коэффициентов $A, B$ путём решения системы
\eqref{abeq} с помощью метода Гаусса\nocite{bahvalov01} (его исходный код приведён на
странице \pageref{gauss.scm-full-listing}):

\input{find-A-B__fundmatrix-solution.scm-tag-listing}

С использованием найденного коэффициента $A$ и приближений из
\procname{build-fundamentals} в функции
\procname{approximate-solution} согласно \eqref{fmsolution}
рассчитывается и приближённое решение исходного уравнения \eqref{deq}
внутри неоднородного слоя:

\input{approximate-solution__fundmatrix-solution.scm-tag-listing}

Функция высшего порядка \procname{iterative-improve} выражает один из
подходов к решению вычислительных задач: улучшать начальное решение
при помощи заданной «функции улучшения» до тех пор, пока его не можно
будет считать «хорошим» в смысле заданной функции оценки решения:

\input{iterative-improve__shared.scm-tag-listing}

Так, в данном методе процедура построения фундаментальных матриц
повторяется до тех пор, пока вычисляемые коэффициенты отражения и
прохождения не будут удовлетворять \eqref{energy}:

\input{get-solution__fundmatrix-solution.scm-tag-listing}

Интересно, что второй метод решения поставленной задачи, о котором
рассказано далее, также определяет сходную функцию
\procname{make-solution} в виде частного случая
\procname{iterative-improve}.

Функцией улучшения является двукратное увеличение числа отрезков
разбиения интервала $[0;a]$ с последующим построением фундаментальных
матриц, решения и нахождения $A$, $B$, а проверка «качества»
полученного решения осуществляется при помощи предиката
\procname{energy-conserves?}:

\input{energy-conserves?__shared.scm-tag-listing}

\clearpage

\subsection{Исходные данные}
\input{statement.scm-initial-data}

В терминах Scheme она задаётся следующим образом:

\input{f__statement.scm-tag-listing}

\subsection{Результаты}

\input{fundmatrix__statement.scm-results}

\includeplot{fundmatrix}{statement.scm}

\clearpage

\subsubsection{Результат вычислений  с альтернативным набором исходных данных}

Метод построения фундаментальных матриц применим в задачах с более
сложными выражениями для $n(x)$ (но требуется, чтобы $n(x)$ была
непрерывной внутри $[0;a]$). Далее представлены результаты работы
с иными, нежели в предыдущем разделе, данными.

\input{exo-statement.scm-initial-data}

\input{fundmatrix__exo-statement.scm-results}

\includeplot{fundmatrix}{exo-statement.scm}

\clearpage

\section{Решение методом последовательных приближений}

\subsection{Описание}

Уравнение \eqref{deq} можно переписать в следующем виде:

\begin{equation}\label{ndeq}
\frac{d^2u}{dx^2} + k^2 u(x) = (k^2 - n(x)) u(x)
\end{equation}

Где $k^2$ — значение функции вне интервала $[0; a]$.

Для такого уравнения известна функция Грина:

\[
G(x, t) = \frac{e^{ik\abs{x-t}}}{2ik}
\]

Так что решение \eqref{ndeq} выписывается в таком виде:

\[
u(x) = \int \limits_{-\infty}^{+\infty} {\frac{e^{ik\abs{x-t}}}{2ik} (k^2-n(t))
  u(t) dt}
\]

В силу свойств функции $n(x)$, при $x < 0$ и $x > a$ имеем $(k^2 -
n(x)) \equiv 0$, так что несобственный интеграл можно заменить на
интеграл в конечных пределах от $0$ до $a$:

\[
u(x) = \int \limits_{0}^{a} {\frac{e^{ik\abs{x-t}}}{2ik} (k^2-n(t))
  u(t) dt}
\]

Кроме того, если раскрыть модуль в показателе $e$ для $x<0$ и $x>a$,
то станет ясно, что выписанное решение содержит волны, уходящие в
$+\infty$ и $-\infty$. Так как постановка задачи содержит волну
$e^{ikx}$, приходящую из минус бесконечности, полное поле во всей
области от $-\infty$ до $+\infty$ имеет следующий вид:

\begin{equation}\label{inteq}
u(x) =  \int \limits_{0}^{a} {\frac{e^{ik\abs{x-t}}}{2ik} (k^2-n(t))
  u(t) dt} + e^{ikx}
\end{equation}

Сокращённо это интегральное уравнение Фредгольма второго рода
записывается в таком виде:

\begin{equation}\label{inteq-short}
u = \op{A}u + f
\end{equation}

Где $\op{A}$ — интегральный оператор $\int_{0}^{a}
{\frac{e^{ik\abs{x-t}}}{2ik} (k^2-n(t)) u(t) dt}$, который действует на
функцию $u(x)$.

При $x < 0$ решение имеет вид:

\[
u(x) = \frac{e^{-ikx}}{2ik} \int \limits_{0}^{a} {e^{ikt} (k^2-n(t))
  u(t) dt} + e^{ikx}
\]

Откуда можно получить выражение для коэффициента $A$:

\begin{equation}\label{int-A}
  A = \frac{1}{2ik} \int \limits_{0}^{a} {e^{ikt} (k^2-n(t)) u(t) dt}
\end{equation}

Аналогично получается следующее выражение и для коэффициента $B$:

\begin{equation}\label{int-B}
  B = \frac{1}{2ik} \int \limits_{0}^{a} {e^{-ikt} (k^2-n(t)) u(t) dt} + 1
\end{equation}

Построение решения $u(x)$ производится при помощи метода
последовательных приближений\nocite{polyanin03}. Решение
\eqref{inteq-short} представим в таком виде:

\[
u = u_0 + u_1 + u_2 + \dotsb
\]

Его подстановка в \eqref{inteq-short} даёт следующее:

\[
 u_0 + u_1 + u_2 + \dotsb = \op{A}u_0 + \op{A}u_1 +
 \op{A}u_2 + \dotsb + f
\]

Если положить $u_0=f, u_1=\op{A}u_0, \dotsc, u_{n+1}=\op{A} u_n$, то
таким выбором последовательных приближений это уравнение будет
тождественно удовлетворено.

После нахождения $u(x)$ коэффициенты $A$ и $B$ вычисляются из
\eqref{int-A} и \eqref{int-B}, соответственно.

\clearpage
\subsection{Реализация}

Для обеспечения работы программы с широким диапазоном исходных данных
— различных функций $n(x)$ — интегрирование \eqref{inteq}
осуществлялось \emph{численно} по формуле Симпсона.

Подинтегральная функция зависит от двух переменных $x$, $t$, а
интегрирование производится по $t$. Используемая функция
\procname{integrate} возвращает $g(x) = \int_a^b{f(x, t) dt}$,
применяя формулу Симпсона так:

\begin{equation}\label{simpson}
  \begin{split}
    \int \limits_a^b {f(x, t) dt} = \frac{h}{3}
    (f(x, a) + 4 (f(x, t_1) + f(x, t_3) + \dotsb + f(x, t_{n-1})) + \\
    + 2 (f(x, t_2) + f(x, t_4) + \dotsb + f(x, t_{n-2})) + f(x, b))
  \end{split}
\end{equation}

Здесь $n = 2k, k \in \mathbb{Z},\ h = (b-a)/n,\ t_k = a+kh$.

\procname{integrate} осуществляет интегрирование функции, переданной
ей в качестве параметра, по второму её аргументу:

\input{integrate__iterative-solution.scm-tag-listing}

Функция \procname{integrate} выражает идею интегрирования функции двух
переменных в общем виде. Для вычисления \eqref{inteq} требуется
использовать приближение функции $u(x)$ и данную по условию функцию
$n(x)$ в выражении $f(x,t) = e^{ik|x-t|}(k^2-n(t))u(t)$, которое и
будет интегрироваться (постоянный множитель $1/2ik$ выносится из под
интеграла). Требуемое сопоставление выполняет функция
\procname{green-transform}:

\input{green-transform__iterative-solution.scm-tag-listing}

Причём эта функция представлена в виде частного случая более общей
процедуры \procname{green-subtransform}, выполняющей преобразование
своих трёх аргументов-функций $u(t), n(t), g(x, t)$ в функцию двух
аргументов $f(x,t) = e^{ik \cdot g(x, t)}(k^2-n(t))u(t)$:

\input{green-subtransform__iterative-solution.scm-tag-listing}

Тогда при $g(x, t) \equiv \abs{x-t}$ получим подинтегральную функцию
из \eqref{inteq}, а при $g(x, t) \equiv -t$ — из выражения для
коэффициента $B$ \eqref{int-B}.

Функция \procname{green-integrate}, используя описанные
\procname{green-transform} и \procname{integrate}, выполняет
интегрирование \eqref{inteq}, принимая функции $u(x), n(x)$ и правый
край отрезка $a$ в качестве аргументов. В сущности,
\procname{green-integrate} представляет в терминах Scheme оператор
$\op{A}$ из \eqref{inteq-short}:

\input{green-integrate__iterative-solution.scm-tag-listing}

В этой функции также осуществляется деление результата интегрирования
на константу $2ik$.

Стоит заметить, что процедура \procname{integrate} (и, следовательно,
\procname{green-integrate}) возвращает не число, а
\emph{функцию}. Вычисление значения этой функции в любой точке $x_0$
порождает свёртку интервала $[0; a]$ в эту точку. Таким образом,
последовательное применение \procname{green-integrate} порождает
функцию, вычисление значения которой породит последовательность
вложенных свёрток на интервале $[0; a]$.

\subsubsection{Применимость метода}
\label{iter-converge}

Сходимость ряда последовательных приближений $u = u_0 + \op{A} u_0 +
\op{A}^2 u_0 \dotsb$ обеспечивается наличием константы $1/2ik$
перед интегралом \eqref{inteq}, если данные $a, n(x), k$ таковы, что
для $M=\max(k^2-n(x_0))$ выполняется:

\begin{equation}\label{iter-converge-cond}
  \frac{Ma}{2k}<1
\end{equation}

В случаях, когда исходные данные не удовлетворяют этому условию, ряд
приближений разойдётся и получить решение не удастся, в то время как
при помощи фундаментальных матриц за приемлемое время решаются
уравнения в системах, где \eqref{iter-converge-cond} не выполняется.

Исходные данные (функция $n(x)$), предлагаемые к использованию в
настоящей курсовой работе, специально подобраны таким образом, чтобы
\eqref{iter-converge-cond} заведомо выполнялись.

\subsubsection{Эффективность реализации}

Высокий уровень применяемой абстракции и простота реализации метода
негативным образом сказываются на эффективности: создаваемая
последовательным действием \procname{green-integrate} композиция
свёрток — функция, вычисление которой в каждой точке требует большого
количества операций.

Сложность \procname{iterative-solve} по элементарным операциям
\procname{integrate} экспоненциально растёт с увеличением числа
применений оператора $\op{A}$ к начальному приближению.

При $u_k(x) = \underbrace{\op{A} \comp \dotsb \comp \op{A}}_k\ \comp\ u_0(x)$ для
вычисления значения $u_k(x_0)$ \emph{в любой точке} $x_0$ потребуется
выполнение $m^k$ операций, где $m$ — постоянное число элементарных
операций сложения, умножения в процедуре \procname{integrate},
зависящее от количества разбиений для интегрирования функции на $[0;
a]$).

Учитывая построение \emph{последовательных} приближений, то есть
последовательное вычисление $A \comp e^{ikx},\ A \comp A \comp
e^{ikx},\ A \comp A \comp A \comp e^{ikx},\dotsc$, количество
выполняемых операций таково, что описываемый подход к реализации
метода очевидно не может реально применяться на практике в расчётных
задачах в силу низкой эффективности.

\subsubsection{Другие подходы к реализации метода}

\paragraph{Символьное интегрирование}

Интеграл \eqref{inteq} может быть вычислен по формуле
Ньютона-Лейбница, если известна первообразная подинтегральной
функции. Её поиск в явном виде составляет задачу символьного
неопределённого интегрирования.

Интегрирование элементарных функций и их простых сочетний является
решаемой алгоритмически задачей, которая, однако, не является
тривиальной (по причинам, в большей степени относящейся к общим
сложностям представления алгебраических выражений на
компьютере). Сложные выражения под знаком интеграла приводят к
непростым заменам, преобразованиям и приёмам. Интеграл может и не
браться в явном виде.

Символьное интегрирование потенциально даёт значительно большие
преимущества, чем численное, но сравнительно сложнее в
реализации. Описанная же ранее численная реализация метода
последовательных приближений представляет собой почти дословный
перевод словесного описания в термины Лиспа.

\paragraph{Использование готового выражения первообразной}

Первообразная $e^{ik\abs{x-t}}(k^2-n(t))u(t)$ по $t$ может быть
получена сторонними средствами — в другой программе или вручную (что
является единственным универсальным и наиболее гибким методом
интегрирования). По затратам на количество операций, выполняемых на
стороне описываемой в настоящей работе программы, этот подход наиболее
эффективен, поскольку сводит задачу нахождения интеграла к совсем
элементарным операциям.

\clearpage
\subsection{Результаты}

Исходные данные к задаче приведены в разделе
\ref{statement.scm-initial-data} на странице
\pageref{statement.scm-initial-data}.

\input{iterative__statement.scm-results}

\includeplot{iterative}{statement.scm}

\clearpage
\section{Решение методом коллокации}
\subsection{Описание}
Рассмотрим данное по условию дифференциальное уравнение \eqref{deq}:

\begin{equation}\label{inteq-eps}
  \epsilon [u(x)] = u'' + n(x) u(x) = 0
\end{equation}

Будем искать его решение в виде

\begin{equation}\label{colloc-approx}
  U_n(x) = \phi_0(x) + \sum_{m=1}^n{C_m \phi_m(x)}
\end{equation}

Здесь $\phi_m(x)$ — система линейно независимых \emph{базисных}
функций. Подставив \eqref{colloc-approx} в \eqref{inteq-eps}, получим
невязку $\epsilon [U_n(x)]$:

\[
\epsilon[U_n(x)]=\left [\phi_0(x)+\sum_{m=1}^n{C_m \phi_m(x)} \right
]''+
n(x){\left [\phi_0(x)+\sum_{m=1}^n{C_m \phi_m(x)} \right ]}
\]

или

\begin{equation}\label{colloc-eps}
\epsilon[U_n(x)]=\psi_0(x)+\sum_{m=1}^n{C_m \psi_m(x)}
\end{equation}

где

\begin{align*}
\psi_0(x)& =\phi_0''(x) + n(x) \phi_0(x) \\
\psi_m(x)& =\phi_m''(x) + n(x) \phi_m(x)
\end{align*}

Метод коллокации предполагает поиск значений $C_m$ таких, что
\eqref{inteq-eps} обращается в 0 в заданной системе из $n$ \emph{точек
  коллокации} $0 \leq x_1 < \dotsb < x_n \leq a$ принадлежащих
рассматриваемому отрезку $[0;a]$:

\begin{equation}
  \epsilon[U_n(x_j)] = 0,\qquad j = 1 \dotso n
\end{equation}

Используя \eqref{inteq-eps} получаем систему из $n$ линейных уравнений
для определения значений $C_m$:

\begin{equation}\label{colloc-system}
  \sum_{m=1}^n{C_m\psi_m(x_j)} = -\psi_0(x_j),\qquad j = 1 \dotso n
\end{equation}

В случае совместности \eqref{colloc-system} найденные коэффициенты
$C_m$ используются в выражении для приближённого решения
\eqref{colloc-approx}.

\subsubsection{О выборе системы базисных функций}
В уравнении \eqref{colloc-eps} члены $\psi_0(x)$ и $\psi_m(x)$
разделены не случайно: при выборе системы $\phi_k$ (а значит и
$\psi_k$) требуется учитывать набор краевых условий \eqref{conds}.

Метод коллокации требует от первой базисной функции $\phi_0(x)$
удовлетворения набору неоднородных краевых условий \eqref{conds}, а от
$\phi_m(x),\ m>0$ — удовлетворения соответствующим однородным краевым
условиями.

\subsection{Реализация}