\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}