\documentclass[12pt,a4paper,oneside,reqno,draft]{amsbook}
\usepackage{amsmath}
\usepackage{amsthm,array,longtable}
\usepackage[T2A]{fontenc}
\usepackage[cp1251]{inputenc}
\usepackage[english,russian]{babel}
\usepackage{amsfonts}
\usepackage{latexsym}


\tolerance 750


\renewcommand{\labelenumi}{\theenumi.}
\renewcommand*{\proofname}{Доказательство}
\DeclareMathOperator*{\grad}{grad}
\DeclareMathOperator*{\Div}{div}
\DeclareMathOperator*{\const}{const}
\newcommand*{\Rd}{\mathbb R^d}
\newcommand*{\sk}{\sum_{|\alpha|=k}}
\newcommand*{\Pc}{\mathcal P}
\newcommand*{\Ac}{\mathcal A}
\newcommand*{\Hc}{\mathcal H}
\newcommand*{\la}{\langle}
\newcommand*{\ra}{\rangle}
\newcommand*{\Sd}{\mathbb S^{d-1}}
\newcommand*{\Bd}{\mathbb B^d}
\newcommand*{\ov}{\overline}


\newtheorem{lemma}{Лемма}[section]
%\newtheorem{theorem}{Теорема}
\newtheorem{theorem}{Теорема}[section]

\newtheorem{proposition}{Предложение}[section]
\newtheorem{corollary}{Следствие}[section]
%\textwidth=160mm
%\usepackage{remreset}
%\makeatletter
%\@removefromreset{section}{chapter}
%\makeatother




\newcounter{exam}[section]
\renewcommand{\theexam}{\arabic{section}.\arabic{exam}}
\newcommand*{\ex}{\par\refstepcounter{exam}%
{\bf Пример \theexam.}\ }
\renewcommand*{\thefigure}{}
\pagestyle{plain}

\begin{document}
\renewcommand{\figurename}{}

\bigskip

\begin{center} Министерство высшего и среднего специального \\
образования РСФСР

\bigskip

Московский авиационный технологический институт\\
им.\ К.Э.~Циолковского

\rule{115mm}{0,4pt}

\vskip10pt

Кафедра ``Высшая математика''

\vskip25pt

\begin{flushright}
Утверждено редакционно-\\
издательским советом института\\
24.10.88
\end{flushright}

\vskip40pt


\large\rm
СПЛАЙН ИНТЕРПОЛЯЦИЯ ФУНКЦИЙ

\vskip30pt

\rm Методические указания к лабораторной работе\\
по курсу ``Высшая математика''

\vskip100pt

\begin{flushright}
Сост.: А.А.~Жуков\\
К.Ю.~Осипенко
\end{flushright}

\vskip 160pt
Москва  1989
\end{center}

\thispagestyle{empty}
\newpage


\setcounter{page}{3}
\section*{ВВЕДЕНИЕ}

Настоящая лабораторная работа посвящена изучению одного из методов аппроксимации функций --- сплайн аппроксимации. Широкое применение сплайны нашли в различных областях техники и технологии в качестве математического аппарата для моделирования поверхностей деталей и агрегатов, сложной формы (обводы различных частей летательных аппаратов, корпусов судов и т.д.). Модели, основанные на сплайн аппроксимации, легли в основу систем автоматизированного проектирования изделий на основе ЭВМ, в которых с помощью аппроксимации сплайнами удалось с достаточной степенью точности решить проблему хранения геометрической информации в числовой форме.

\bigskip

\begin{center}
Порядок выполнения работы
\end{center}

1. Ознакомиться с методом сплайн интерполяции.

2а. (В случае выполнения работы на ПЭКВМ). Для своего варианта составить вспомогательную таблицу 2. Вычислить значение сплайна в заданной точке.

2б. (В случае выполнения работы на ЭВМ). Ознакомиться со стандартной программой для вычисления значений кубического сплайна на равномерной сетке. Для своего варианта составить подпрограмму обращения к стандартной и посчитать значение сплайна в заданной точке.

3. Составить отчет о работе.

4. Защита работы.

\section{ОПИСАНИЕ МЕТОДА СПЛАЙН ИНТЕРПОЛЯЦИИ}

Пусть отрезок $[a,b]$ разбит на $N$ равных отрезков $[x_i,x_{i+1}]$, где $x_i=a+ih$, $i=0,1,\ldots,N-1$; $x_N=b$; $h=(b-a)/N$. Сплайном называется функция, которая вместе с несколькими производными непрерывна на всем заданном отрезке $[a,b]$, а на каждом частичном отрезке $[x_i,x_{i+1}]$ является многочленом фиксированной степени.

На практике широкое распространение получили кубические сплайны $S_3(x)$, являющиеся непрерывными на $[a,b]$ функциями вместе со своими первой и второй производными, совпадающие с многочленом третьей степени на каждом отрезке $[x_i,x_{i+1}]$, $i=0,1,\ldots,N-1$. Кубический сплайн, принимающий в узлах $x_i$ те же значения $y_i$, что и некоторая функция $f$, называется интерполяционным. Он применяется для аппроксимации функции $f$ на отрезке $[a,b]$ и ее производных до третьего порядка.

Для построения интерполяционного кубического сплайна положим $M_i=S_3''(x_i)$, $i=0,1,\ldots,N$ и заметим, что в силу линейности функции $S_3''(x)$ на каждом из отрезков $[x_i,x_{i+1}]$
\begin{equation}\label{1}
S_3''(x)=M_i\frac{x_{i+1}-x}h+M_{i+1}\frac{x-x_i}h,\quad x\in[x_i,x_{i+1}].
\end{equation}
Интегрируя дважды обе части равенства \eqref{1} и вычисляя константы из условий $S_3(x_i)=f(x_i)$, $S_3(x_{i+1})=f(x_{i+1})$, получим
\begin{multline}\label{2}
S_3(x)=M_i\frac{(x_{i+1}-x)^3}{6h}+M_{i+1}\frac{(x-x_i)^3}{6h}+\left(y_i-\frac16M_ih^2\right)
\frac{x_{i+1}-x}h\\
+\left(y_{i+1}-\frac16M_{i+1}h^2\right)\frac{x-x_i}h,\quad x\in[x_i,x_{i+1}].
\end{multline}
Из равенства \eqref{2} найдем односторонние производные
\begin{equation}\label{3}
\begin{gathered}
S_3'(x_i-0)=\frac h6M_{i-1}+\frac h3M_i+\frac{y_i-y_{i-1}}h,\\
S_3'(x_i+0)=-\frac h3M_i-\frac h6M_{i+1}+\frac{y_{i+1}-y_i}h.
\end{gathered}
\end{equation}
Поскольку функция $S_3'(x)$ непрерывна на $[a,b]$, то эти производные должны быть равными. Таким образом, получается система уравнений
\begin{equation}\label{4}
\frac16M_{i-1}+\frac{2h}3M_i+\frac h6M_{i+1}=\frac{y_{i-1}-2y_i+y_{i+1}}h,\quad
i=1,2,\ldots,N-1.
\end{equation}
В этой системе $N-1$ уравнение. Чтобы однозначно определить значения $M_0,M_1,\ldots,M_N$, надо задать два дополнительных условия, так называемые ``краевые условия''.

Краевые условия могут задаваться различным образом. Рассмотрим следующие случаи.

1. На концах известны значения вторых производных $M_0=y_0''$, $M_N=y_N''$, в частности, $M_0=M_N=0$.

2. На концах известны наклоны $S_3'(a)=y_0'$, $S_3'(b)=y_N'$. Тогда из равенств \eqref{3} получаем
\begin{gather*}
\frac h3M_0+\frac h6M_1=\frac{y_1-y_0}h-y_0',\\
\frac h6M_{N-1}+\frac h3M_N=y_N'-\frac{y_N-y_{N-1}}h.
\end{gather*}

Оба этих случая могут быть объединены одним видом краевых условий
\begin{equation}\label{5}
2M_0+\lambda M_1=d_0,\quad\lambda M_{N-1}+2M_N=d_N,
\end{equation}
где в случае 1) $\lambda=0$, $d_0=2y_0''$, $d_N=2y_N''$, а в случае 2) $\lambda=1$, $d_0=6((y_1-y_0)/h-y_0')/h$, $d_N=6(y_N'-(y_N-y_{N-1})/h)/h$.

Тем самым система \eqref{4} с учетом краевых условий \eqref{5} принимает вид
\begin{equation}\label{6}
\begin{pmatrix}
2&\lambda&0&\ldots&0&0&0\\
1&4&1&\ldots&0&0&0\\
0&1&4&\ldots&0&0&0\\
\hdotsfor{7}\\
0&0&0&\ldots&4&1&0\\
0&0&0&\ldots&1&4&1\\
0&0&0&\ldots&0&\lambda&2
\end{pmatrix}
\begin{pmatrix}
M_0\\M_1\\M_2\\\vdots\\M_{N-2}\\M_{N-1}\\M_N
\end{pmatrix}=
\begin{pmatrix}
d_0\\d_1\\d_2\\\vdots\\d_{N-2}\\d_{N-1}\\d_N
\end{pmatrix},
\end{equation}
где
$$d_i=\frac6{h^2}(y_{i+1}-2y_i+y_{i-1}),\quad i=1,2,\ldots,N-1.$$

Для решения системы \eqref{6} образуем вспомогательные величины
\begin{gather*}
q_0=-\lambda/2,\quad p_N=\lambda q_{N-1}+2,\\
p_i=q_{i-1}+4,\quad q_i=-1/p_i,\quad i=1,2,\ldots,N-1,\\
u_0=d_0/2,\quad u_i=(d_i-u_{i-1})/p_i,\quad i=1,2,\ldots,N-1,\\
u_N=(d_N-\lambda u_{N-1})/p_N.
\end{gather*}
Далее последовательно определяем $M_N,M_{N-1},\ldots,M_0$
$$M_N=u_N,\quad M_i=q_iM_{i+1}+u_i,\quad i=N-1,\ldots,0.$$
Зная величины $M_0,\ldots,M_N$, по формулам \eqref{2} можно вычислить значение сплайна $S_3(x)$ в любой точке.

\section{СХЕМА ЗАПИСИ ВЫЧИСЛЕНИЙ}

Величины $p_i$ и $q_i$ не зависят от $y_i$, а зависят лишь от $\lambda$. Поэтому они могут быть вычислены заранее для $\lambda=0,1$. Начиная с $N=6$, эти величины с точностью до 5-го знака не изменяются. При данном $N$ надо взять значения $p_i$ и $q_i$, $i=0,\ldots,N-1$ из таблицы~1 для заданного $\lambda=0$ или $1$, а значение $p_N$ из столбца поправок для $\lambda=1$ и $p_N=2$ для $\lambda=0$.

\centerline{Таблица~1}

\bigskip

\begin{center}
\begin{tabular}{|l|l|l|l|l|l|l|}
\hline
&\multicolumn{2}{|c|}{$\lambda=0$}&\multicolumn{2}{|c|}{$\lambda=1$}&
\multicolumn{2}{|c|}{$\lambda=1$}\\
\hline
\multicolumn{1}{|c|}{$i$}&\multicolumn{1}{|c|}{$p_i$}&\multicolumn{1}{|c|}{$q_i$}
&\multicolumn{1}{|c|}{$p_i$}&\multicolumn{1}{|c|}{$q_i$}&
\multicolumn{1}{|c|}{\raisebox{-1,5pt}{$N$}}&
\multicolumn{1}{|c|}{$p_N$}\\
\hline
0&--&\hphantom-0&--&\hphantom-0,5&&\\
1&4&-0,25&3,5&-0,28571&1&1,6\\
2&3,75&-0,26667&3,71429&-0,26923&2&1,71429\\
3&3,73333&-0,26788&3,73077&-0,26804&3&1,73077\\
4&3,73214&-0,26794&3,73196&-0,26796&4&1,73196\\
5&3,73206&-0,26795&3,73204&-0,26795&5&1,73204\\
6&3,73206&-0,26795&3,73205&-0,26795&6&1,73205\\
7&3,73206&-0,26795&3,73205&-0,26795&7&1,73205\\
\hline
\end{tabular}
\end{center}

Вычисления удобно расположить в следующей таблице.

\bigskip

\centerline{Таблица~2}

$$\begin{array}{|c|c|c|c|c|c|c|}
\hline
i&p_i&q_i&y_i&d_i&u_i&M_i\\
\hline
0&$--$&q_0&y_0&d_0&u_0&M_0\\
1&p_1&q_1&y_1&d_1&u_1&M_1\\
2&p_2&q_2&y_2&d_2&u_2&M_2\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\
N-1&p_{N-1}&q_{N-1}&y_{N-1}&d_{N-1}&u_{N-1}&M_{N-1}\\
N&p_N&q_N&y_N&d_N&u_N&M_N\\
\hline
\end{array}$$

\bigskip

Пример.

Построить кубический интерполяционный сплайн для функции $y=\sin x$ на отрезке $[0,\pi/2]$ для $N=5$ с краевыми условиями $y_0''=0$, $y_5''=-1$ ($\lambda=0$). Вычислить значение сплайна в точке $x=\pi/4$ и сравнить с точным значением функции.

Составим таблицу, аналогичную таблице~2.

$$\begin{array}{|c|l|l|l|l|l|l|}
\hline
i&\multicolumn{1}{c|}{p_i}&\multicolumn{1}{c|}{q_i}&\multicolumn{1}{c|}{y_i}&
\multicolumn{1}{c|}{\alpha_i}&\multicolumn{1}{|c|}{u_i}&\multicolumn{1}{|c|}{M_i}\\
\hline
0&\multicolumn{1}{c|}{$--$}&\hphantom{-}0&0&\hphantom{-}0&\hphantom{-}0&\hphantom{-}0\\
1&4&-0,25&0,30902&-1,83898&-0,45975&-0,31154\\
2&3,75&-0,26667&0,58779&-3,49801&-0,81020&-0,59286\\
3&3,73333&-0,26786&0,80902&-4,81417&-1,07249&-0,81502\\
4&3,73214&-0,26794&0,95106&-5,65980&-1,22914&-0,96120\\
5&2&\multicolumn{1}{c|}{$--$}&1&-2&-1&-1\\
\hline
\end{array}$$

\bigskip

Для вычисления сплайна в точке $x=\pi/4$ воспользуемся формулой~\eqref{2} при $i=2$.
\begin{multline*}
S_3\left(\frac\pi4\right)=-0,59286\frac{\left(\dfrac{3\pi}{10}-\dfrac\pi4\right)^3}
{\dfrac{6\pi}{10}}-0,81502\frac{\left(\dfrac\pi2-\dfrac{2\pi}{10}\right)^3}
{\dfrac{6\pi}{10}}
+\biggl(0,58779\\+\frac{0,59286}6\left(\frac\pi{10}\right)^2\biggr)\frac{\dfrac{3\pi}{10}-\dfrac\pi4}
{\dfrac\pi{10}}+\left(0,80902+\frac{0,81502}6\left(\frac\pi{10}\right)^2\right)
\frac{\dfrac\pi4-\dfrac{2\pi}{10}}{\dfrac\pi{10}}\\
=-\left(\frac\pi{10}\right)^2\frac{0,59286+0,81502}{48}+\frac{0,58779+0,80902}2\\
+\left(\frac\pi{10}\right)^2\frac{0,59286+0,81502}{12}=
\left(\frac\pi{10}\right)^2(-0,02933+0,11732)+0,69841\\
=0,70709.
\end{multline*}

Точное значение $y\left(\dfrac\pi4\right)=\dfrac{\sqrt2}2=0,70711$.

\section{СТАНДАРТНАЯ ПОДПРОГРАММА ДЛЯ ВЫЧИСЛЕНИЯ ЗНАЧЕНИЙ КУБИЧЕСКОГО СПЛАЙНА}

Для вычисления значений интерполяционного кубического\break сплайна по равномерной сетке можно использовать стандартную подпрограмму $SPL3$. Обращение к этой подпрограмме имеет вид
$$CALL\ SPL3(A,B,DA,DB,M,Y,RL,X,YS,K,Q,U).$$

Формальные параметры имеют следующий смысл:

\bigskip

\begin{tabular}{c@{\ ---}lp{9,2cm}}
$A$&&левый конец отрезка;\\
$B$&&правый конец отрезка;\\
$DA$, $DB$&&соответственно, $d_0$ и $d_N$ (см.\ \eqref{5});\\
$M$&&число точек равномерного разбиения отрезка $[A,B]$ (оно равно $N+1$);\\
$Y$&&массив значений приближаемой функции в равномерной сетке отрезка $[A,B]$ (размерности $M$);\\
$RL$&&параметр $\lambda$ (см.\ \eqref{5});\\
$X$&&массив значений аргумента, в которых требуется посчитать значения сплайна (размерности $K$);\\
$YS$&&массив значений сплайна в точках массива $X$ (размерности $K$);\\
$X$&&размерность массивов $X$ и $YS$;\\
$Q$, $U$&&рабочие массивы размерности $M$.
\end{tabular}

\bigskip

Подпрограмма $SPL3$ имеет следующий вид

\newcolumntype{L}{>{$}l<{$}}

\begin{longtable}{r|c|L}
&&SUBROUTINE\ SPL3(A,B,DA,DB,M,Y,RL,X,\\
&*&YS,K,Q,U)\\
&&DIMENSION\ Y(M),X(K),YS(K),Q(M),U(M)\\
&&M1=M-1\\
&&H=(B-A)/M1\\
&&H1=6./H**2\\
&&Q(1)=-.5*RL\\
&&U(1)=-.5*DA\\
&&DO\ 1\ I=2,M1\\
&&Q(I)=-1./(Q(I-1)+4.)\\
&&D=H1*(Y(I+1)-2*Y(I)+Y(I-1))\\
1&&U(I)=(D-U(I-1))/(Q(I-1)+4.)\\
&&U(M)=(DB-RL*U(M1))/(RL*Q(M1)+2.)\\
&&DO\ 2\ I=1,M1\\
&&L=M-I\\
2&&U(L)=Q(L)*U(L+1)+U(L)\\
&&DO\ 3\ I=1,K\\
&&L=M1\\
&&DO\ 4\ J=2,M1\\
&&IF(X(I).GT.A+(J-1)*H)\ GOTO\ 4\\
&&L=J-1\\
&&GOTO\ 5\\
4&&CONTINUE\\
5&&X2=A+L*H\\
&&X1=X2-H\\
3&&YS(I)=(U(L)*(X2-X(I))**3+U(L+1)*(X(I)-\\
&*&X1)**3/6./H+(Y(L)-U(L)/H1)*(X2-X(I))\\
&*&/H+(Y(L+1)-U(L+1)/H1)*(X(I)-X1)/H\\
&&RETURN\\
&&END
\end{longtable}

\section{ВАРИАНТЫ ИНДИВИДУАЛЬНЫХ ЗАДАНИЙ}

В каждом варианте для заданной функции $f(x)$ построить на указанном отрезке $[a,b]$ интерполяционный кубический сплайн для $N=5$ с краевыми условиями:
$$y_0''=f''(a),\quad y_5''=f''(b).$$

Вычислить значение сплайна в середине отрезка $\dfrac{a+b}2$ и сравнить его с точным значением функции.

\newcolumntype{M}{>{$}l<{,$}}
\newcolumntype{N}{>{$[}l<{]$}}

\begin{longtable}{rMN}
1.&\tg x^2&0,\pi/3\\
2.&\arctg\sqrt x&1/2,1\\
3.&\ln(x^2+1)&-2,-1\\
4.&e^{x-1}&0,1\\
5.&\arctg x&0,1\\
6.&\ctg(x+1)&-1/2,0\\
7.&\ln(sin^2x)&\pi/3,\pi/2\\
8.&\cos(x^2-1)&0,\pi/2\\
9.&2^{\sin x}&0,\pi/2\\
10.&\arctg\sqrt{x-1}&2,3\\
11.&\ln(x^2+x)&1/2,1\\
12.&1/\ln x&4,5\\
13.&\arctg(x^2-1)&0,1\\
14.&1/\arcsin x&1/2,1\\
15.&\sin\sqrt x&1/2,1\\
16.&\tg(1/x^2)&2,3\\
17.&1/\sqrt{x^2+4}&0,1\\
18.&\arcsin x^2&0,1/2\\
19.&e^{x^2}&1,3/2\\
20.&e^x/\sqrt x&1,2\\
21.&\sqrt{x^2+x-1}&2,3\\
22.&\ln(\arctg x)&1,2\\
23.&1/((x^2+1)&-1/2,1/2\\
24.&3^{x^2+1}&-1,0\\
25.&1/\arctg x&1,2\\
26.&\cos(1/x)&2,3\\
27.&\ln(\arcsin x)&1/2,1\\
28.&(2x+1)/x^3&1,2\\
29.&e^{1/x}&1,2\\
30.&\sqrt{3-2x+x^2}&-3,-2.5
\end{longtable}

\bigskip

\centerline{Основная литература}
\begin{itemize}
\item[1.] Бахвалов~Н.С., Жидков~Н.П., Кобельков~Г.М. Численные методы. --- М.: Наука, 1987.
\item[2.] Калиткин~Н.Н., Численные методы. --- М.: Наука, 1978.
\end{itemize}

\bigskip

\centerline{Дополнительная литература}
\begin{itemize}
\item[1.] Альберг~Дж., Нильсон~Э., Уолш~Дж. Теория сплайнов и ее приложение. --- М.: Мир, 1972.
\item[2.] Завьялов~Ю.С., Квасов~Б.И., Мирошниченко~В.Л. Методы сплайн-функций. --- М.: Наука, 1980.
\item[3.] Форсайт~Дж., Мальком~М., Моулер~К. Машинные методы математических вычислений. --- М.: Мир, 1980.
\end{itemize}
\end{document}





