\documentclass[12pt,draft,a4paper]{amsart}
\usepackage{amsmath,amsthm,array,longtable}
\usepackage[T2A]{fontenc}
\usepackage[cp1251]{inputenc}
\usepackage[english,russian]{babel}
\usepackage{amsfonts}
\usepackage{latexsym}
\usepackage{cite}
%\usepackage{abstract}
%\usepackage{srctex}
\tolerance 1000

\makeatletter
\renewcommand{\@biblabel}[1]{#1.}
\makeatother

\newtheorem{theorem}{Теорема}
\newtheorem*{lemma}{Лемма}
\newtheorem{proposition}{Предложение}
\newtheorem*{corollary}{Следствие}
\renewcommand{\thesection}{\S \arabic{section}}
\newtheorem*{theorem1}{Теорема~$\mathbf1'$}

\DeclareMathOperator*{\infp}{inf\vphantom p}
\DeclareMathOperator*{\vraisup}{vraisup}
\DeclareMathOperator*{\IM}{Im}
\DeclareMathOperator*{\RE}{Re}
\DeclareMathOperator*{\gr}{gr}
\DeclareMathOperator{\sn}{sn}
\DeclareMathOperator*{\spann}{span}
\DeclareMathOperator{\cn}{cn}
\DeclareMathOperator{\ctn}{ctn}
\DeclareMathOperator{\dn}{dn}
\DeclareMathOperator{\sign}{sign}
\DeclareMathOperator{\arth}{arth}
\DeclareMathOperator{\thh}{th}
\DeclareMathOperator{\Aff}{Aff}
\DeclareMathOperator{\co}{co}
\DeclareMathOperator{\bco}{bco}
\DeclareMathOperator*{\ulim}{\underline{\lim}}

\newcommand*{\ov}{\overline}
\newcommand*{\ei}{e^{i\theta}}
\newcommand*{\la}{\langle}
\newcommand*{\ra}{\rangle}
\newcommand*{\wt}{\widetilde\Theta}
\newcommand*{\wl}{\widetilde l}
\newcommand*{\pph}[1]{\Phi_{\lambda,#1,\beta}}
\newcommand*{\pr}{\texttt{\textvisiblespace}}
\newcommand*{\wy}{\widehat y}

\begin{document}
\begin{center}
К.~Ю.~Осипенко\\
\bigskip

ПОДПРОГРАММА ВЫЧИСЛЕНИЯ ПЕРВОЙ И ВТОРОЙ ПРОИЗВОДНОЙ ФУНКЦИИ ОДНОГО ПЕРЕМЕННОГО МЕТОДОМ СГЛАЖИВАЮЩИХ СПЛАЙНОВ
\end{center}

\bigskip

Подпрограмма $SM${\it\O\O\/}$TH$ предназначена для вычисления значений функции действительного переменного вместе с ее первой и второй производной по приближенным значениям функции в узлах произвольной сетки. В подпрограмме используется метод сглаживающих сплайнов \cite{1}, \cite{2}.

Подпрограмма написана на языке ФОРТРАН \cite{3}.

\refstepcounter{section}
\section*[\ ]{\S\arabic{section}. Обращение к подпрограмме}

Обращение к подпрограмме имеет вид:
$$CALL\ SM\mbox{\it\O\O\/}TH(X,Y,P,D,E,K,N,D1,R\mbox{\it\O},Y0,Y1,Y2,R,A)$$
Здесь формальные параметры имеют следующий смысл:
\begin{itemize}
\item[$X$ ---] вещественный одномерный массив (размерности $N$) значений аргумента функции.
\item[$Y$ ---] вещественный одномерный массив (размерности $N$) значений функции в узлах сетки.
\item[$P$ ---] вещественный одномерный массив (размерности $N$) значений весовых коэффициентов.
\item[$D$ ---] числовой параметр, определяющий начало итерационного процесса (см.\ описание метода в \S~2).
\item[$E$ ---] числовой параметр, характеризующий среднеквадратическую точность задаваемых значений функции в узлах сетки.
\item[$K$ ---] условное число, принимающее целые значения $0$ или $1$ и задаваемое пользователем.
\end{itemize}

Если $K=1$, то в процессе работы подпрограммы $E$ присваивается значение $E*R\mbox{\it\O}$, где $R\mbox{\it\O}$ --- среднеквадратичное отклонение от прямой, получаемой по методу наименьших квадратов.

\begin{itemize}
\item[$N$ ---] число точек сетки.
\item[$D1$ ---] значение параметра $D$ при выходе из итерационного процесса (см.\ описание метода в \S~2).
\item[$R\mbox{\it\O}$ ---] значение достигнутой точности приближения.
\item[$Y0$ ---] вещественный одномерный массив (размерности $N$) значений сглаживающего сплайна в узлах сетки.
\item[$Y1$ ---] вещественный одномерный массив (размерности $N$) значений первой производной сглаживающего сплайна в узлах сетки.
\item[$Y2$ ---] вещественный одномерный массив (размерности $N$) значений второй производной сглаживающего сплайна в узлах сетки.
\item[$R$ ---] значение сглаживающего функционала
$$R=\int_{x_1}^{x_N}\left[y''(x)\right]^2\,dx.$$
\item[$A$ ---] вещественный двумерный массив (размерности $N\times7$) рабочих ячеек.
\end{itemize}

\refstepcounter{section}
\section*[\ ]{\S\arabic{section}. Описание метода построения сглаживающего сплайна}

Ставится задача нахождения функции (сглаживающего сплайна), которая минимизирует функционал
$$R=\int_{x_1}^{x_N}\left[y''(x)\right]^2\,dx$$
на классе функций, имеющих непрерывную вторую производную на отрезке $[x_1,x_N]$ и удовлетворяющую условию
$$\sum_{i=1}^Np_i(y(x_i)-y_i)^2\le E^2,$$
где $y_i$ --- задаваемые значения табличной функции $y(x)$ в узлах сетки $x_i$ ($i=1,2,\ldots,N$), $p_i>0$ --- весовые коэффициенты, $E$ --- абсолютная или относительная точность задаваемых значений функции $y_i$, определяемая параметром $K$: при $K=0$ величина $E$ определяет абсолютную точность, при $K=1$ величина $E$, $0<E<1$, определяет среднеквадратичную точность $E:=E\rho(0)$, где $\rho(0)$ --- среднеквадратическая погрешность приближения заданных значений прямой линией, определяемая программно.

Введем обозначения: $h_i=x_{i+1}-x_i$, $i=1,\ldots,N-1$, и определим матрицы:

\begin{gather*}
Q^T=\arraycolsep=0.06em\begin{pmatrix}
\frac1{h_1}&-\left(\frac1{h_1}+\frac1{h_2}\right)&\frac1{h_2}&0&\ldots&0&0&0\\
0&\frac1{h_2}&-\left(\frac1{h_2}+\frac1{h_3}\right)&\frac1{h_3}&\ldots&0&0&0\\
\hdotsfor{8}\\
0&0&0&0&\ldots&\frac1{h_{N-2}}&-\left(\frac1{h_{N-2}}+\frac1{h_{N-1}}\right)&
\frac1{h_{N-1}}
\end{pmatrix}\\
T=\begin{pmatrix}\dfrac13(h_1+h_2)&\dfrac{h_2}6&0&\ldots&0&0\\
\dfrac{h_2}6&\dfrac13(h_2+h_3)&\dfrac{h_3}6&\ldots&0&0\\
\hdotsfor{6}\\
0&0&0&\ldots&\dfrac{h_{N-2}}6&\dfrac13(h_{N-2}+h_{N-1})
\end{pmatrix}
\end{gather*}
$$P=\begin{pmatrix}
p_1&&&\\
&p_2&\raisebox{7pt}[0pt][0pt]{\hbox to 0pt{\hphantom{0}$0$\hss}}&\\
&\raisebox{-7pt}[0pt][0pt]{\hbox to 0pt{\hss$0$\hphantom{0}}}&\ddots&\\
&&&p_N
\end{pmatrix}.$$


Положим
$$\|z\|_p^2=\sum_{i=1}^Np_iz_i^2$$
и
$$\rho(d)=\|P^{-1}Q(Q^TP^{-1}Q+dT)^{-1}Q^T\wy\|_p,$$
где $\wy=(y_1,y_2,\ldots,y_N)^T$, $d>0$ --- вещественный параметр.

В работах \cite{1}, \cite{2} показано, что решение поставленной ранее экстремальной задачи является кубический сплайн $y_0(x)$, значения которого в узлах сетки $\{x_i\}$ определяются компонентами вектора $\wy_d$:
$$\wy_d=\wy-P^{-1}Q\ov s,\quad(Q^TP^{-1}Q+dT)\ov s=Q^T\wy,$$
где вектор $\ov s=(s_2,s_3,\ldots,s_{N-2})$, а параметр $d$ является (единственным) решением скалярного уравнения:
$$\rho(d)=E.$$
Для решения этого уравнения применяется итерационный метод касательных Ньютона, при этом начальное значение параметра $d=D$, а вычисленное программой --- $D1$. Всегда можно положить $D=0$, однако, если проводится серия расчетов при близких значениях $E$, то для ускорения времени счета целесообразно задавать значение $D=D1$, полученному на предыдущем шаге.

Значение первой $y_1(x_i)$ и второй $y_2(x_i)$ производной сглаживающего сплайна в узловых точках сеток определяются по формулам:
\begin{enumerate}
\item[I)]$y_1(x_i)=\dfrac{y_0(x_{i+1})-y_0(x_i)}{h_i}-\dfrac{2y_2(x_i)+y_2(x_{i+1})}6h_i$,
\item[]\hspace{229pt}$i=1,2,\ldots,N-1$,
\item[]$y_1(x_N)=y_1(x_{N-1})+\dfrac{y_2(x_{N-1})}2h_{N-1}$;
\item[II)]$y_2(x_i)=ds_i$, $i=2,\ldots,N-1$, $y_2(x_1)=y_2(x_N)=0$.
\end{enumerate}
Значение $R=\displaystyle\int_{x_1}^{x_N}\left[y_2(x)\right]^2\,dx$ определяется равенством
$$R=(T\ov s,\ov s).$$

\refstepcounter{section}
\section*[\ ]{\S\arabic{section}. Тестовый пример}

Пусть $x_i=\dfrac{i-1}{10}$, $p_i=1$, $i=1,\ldots,30$, а $y_i$ --- значения $\sin x_i$ с точностью до трех знаков \rule{0pt}{10pt}

$$\arraycolsep=2,84em\begin{array}{rlc}
y_1=.000&y_{11}=.841&y_{21}=.909\\
y_2=.100&y_{12}=.891&y_{22}=.863\\
y_3=.199&y_{13}=.932&y_{23}=.808\\
y_4=.296&y_{14}=.964&y_{24}=.746\\
y_5=.389&y_{15}=.985&y_{25}=.675\\
y_6=.479&y_{16}=.997&y_{26}=.598\\
y_7=.565&y_{17}=1.000&y_{27}=.516\\
y_8=.644&y_{18}=.992&y_{28}=.427\\
y_9=.717&y_{19}=.974&y_{29}=.335\\
y_{10}=.783&y_{20}=.946&y_{30}=.239\\
\end{array}$$
Было задано $E=\sqrt{2.5}\times10^{-3}$, $D=0$. Обращение к подпрограмме имело вид:
$$CALL\ SM\mbox{\it\O\O\/}TH(X,Y,P,0.,E,0,30,D1,R\mbox{\it\O},Y0,Y1,Y2,R,A)$$

Были получены следующие результаты (приводятся с точностью до пяти знаков после запятой).

\tabcolsep=2.36em\begin{longtable}{rrrr}
\multicolumn{1}{c}{$x_i$}&\multicolumn{1}{c}{Y0}&\multicolumn{1}{c}{Y1}&
\multicolumn{1}{c}{Y2}\\
\hline
0.0&0.00030&0.99965&0.00000\\
0.1&0.10011&0.99513&-0.09043\\
0.2&0.19897&0.97985&-0.21524\\
0.3&0.29568&0.95259&-0.32996\\
0.4&0.38926&0.91863&-0.34921\\
0.5&0.47921&0.87881&-0.44723\\
0.6&0.56459&0.82595&-0.60993\\
0.7&0.64407&0.76305&-0.64797\\
0.8&0.71704&0.69531&-0.70686\\
0.9&0.78292&0.62113&-0.77669\\
1.0&0.84107&0.54123&-0.82134\\
1.1&0.89098&0.45578&-0.88777\\
1.2&0.93202&0.36401&-0.94750\\
1.3&0.96357&0.26602&-1.01240\\
1.4&0.98522&0.16799&-0.94822\\
1.5&0.99727&0.07306&-0.95032\\
1.6&0.99969&-0.03621&-1.03504\\
1.7&0.99191&-0.12921&-1.02508\\
1.8&0.97392&-0.22982&-0.98710\\
1.9&0.94611&-0.32548&-0.92605\\
2.0&0.90898&-0.41666&-0.89750\\
2.1&0.86288&-0.50463&-0.86194\\
2.2&0.80823&-0.58727&-0.79091\\
2.3&0.74555&-0.66626&-0.78882\\
2.4&0.67521&-0.73824&-0.65071\\
2.5&0.59826&-0.79957&-0.57597\\
2.6&0.51542&-0.85731&-0.57890\\
2.7&0.42708&-0.90652&-0.40520\\
2.8&0.33465&-0.93955&-0.25534\\
2.9&0.23985&-0.95231&-0.00000\\
\hline
\end{longtable}
$D1=3020,98091$, $R\mbox{\it\O}=0.00158$, $R=1.55938$.

\bigskip

\renewcommand{\refname}{ЛИТЕРАТУРА}
\begin{thebibliography}{11}

\bibitem{1} Морозов~В.А. О задаче дифференцирования и некоторых алгоритмах приближения экспериментальной информации. Сб. работ ``Вычислительные методы и программ.'', вып.~XIV, изд. Моск. ун-та, М., 1970, 46--62.
\selectlanguage{english}
\bibitem{2} Reinsch~C.H. Smoothing by splines functions. Numer. Math., 1967, 10,\selectlanguage{russian} \No~3, 177--183.
\bibitem{3} Язык ФОРТРАН. Под ред. В.П.~Ширикова, ВЦ МГУ, М., 1970.

\end{thebibliography}
\end{document}