Сплайн-интерполяция

ПРИДНЕСТРОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

 им. Т.Г.  Шевченко

 

Рыбницкий филиал

 

кафедра физики, математики и информатики 

 

 

  

 

 

 

 

 

 

КУРСОВАЯ РАБОТА

 

по дисциплине «Вычислительная математика»

 

Тема: «Сплайн-интерполяция» 

 

 

 

 

Выполнил:

студент  II курса 220 гр,

специальности

«ПОВТ и АС»

 

 

Проверил:

ст .преподаватель

Балан  Л. А.  

 

 

 

 

 

 

 

Рыбница

2011 г.

 

 

Оглавление

Введение 2

Глава  I. Теоретические основы интерполирования кубическим сплайном. 4

1.1 Постановка задачи интерполирования кубическими сплайнами. 4

1.2 Алгоритм построения интерполяционного кубического сплайна. 5

1.3 Обзор методов решения СЛАУ с трехдиагональными матрицами методом фронтальной прогонки. 9

1.4 Численный пример. 12

Глава II. Практическая реализация интерполирования кубическим сплайном 19

2.1 Реализация в С++. 19

2.2 Реализация в Exel. 20

2.3 Реализация в MathCad. 22

2.4 Сравнение  полученных результатов. 23

2.5 Блок-схема  программной реализации сплайн  интерполяции. 26

Заключение. 27

Список  литературы 28

Приложение. 29

 

 

 

 

 

 

 

 

 

Введение

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

Во всех этих случаях используются методы приближенного, в первую очередь  численного решения. Методы численного решения математических задач всегда составляли неотъемлемую часть математики и неизменно входили в содержание естественно-математического и инженерного  образования.

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

Технологическая цепочка вычислительного  эксперимента включает в себя следующие  этапы:

• построение математической модели исследуемого объекта (сюда же относится и анализ модели, выяснение корректности поставленной математической задачи;
• построение вычислительного алгоритма - метода приближенного решения поставленной задачи и его обоснование;
• программирование алгоритма на ЭВМ и его тестирование;
• проведение серии расчетов с варьированием определяющих параметров исходной задачи и алгоритма;
• анализ полученных результатов;

Целью данной курсовой работы является разработка программного продукта, реализующего решение  интерполяционного кубического  сплайна методом фронтальной  прогонки, а так же методом Гаусса.

Для достижения цели были поставлены  следующие задачи:

1. Изучить теоретический материал по теме: «Интерполяционные кубические сплайны».
2. Применить изученный материал к разработке программного продукта.
3. Написать программу на языке С++ для решения данного кубического сплайна методом фронтальной рпогонки.
4. Проверить решение с помощью приложения MathCad, Excel.
5. Сравнить результаты, полученные в трех пакетах.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Глава  I. Теоретические основы интерполирования кубическим сплайном.

1. Постановка задачи интерполирования кубическими сплайнами.

Пусть на отрезке [a,b] задана сеточная функция у₁ = f(x_i) на сетке Ω_п. Требуется восполнить ее функцией S_m(x), где т — степень многочлена, кусочно-глобальным способом.

Дадим сначала упрощенное определение  сплайна, которое затем уточним.

Сплайн-функцией или сплайном называется совокупность S_m,i(x)- алгебраических многочленов степени m (звеньев), определенных на частичных отрезках [x_i,x_i+1] , i= , и соединенных вместе по всем частичным отрезкам так, чтобы можно было составить многозвенную функцию

S_m(x)=    определенную и непрерывную на всем отрезке [a,b] вместе со всеми своими производными до некоторого их порядка p=1,2,…..

Разность между р и наибольшим порядком производной, непрерывной на отрезке [а, b], определяет дефект  сплайна q.

Условиями согласования звеньев S_m,i(x) сплайна с исходной функцией                       y_i=f(x_i) на соответствующем частичном отрезке [х_i,х_i+1] называются условия, накладываемые на невязки дифференциального и интегрального типов   , и использующиеся для вывода формулы одного звена сплайна на указанном частичном отрезке.

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

 

 

 

2. Алгоритм построения интерполяционного кубического сплайна.

Первый этап. Для алгебраического многочлена (полинома)

S_3,i(x) =  a_0,i + a_1,i(x-x_i) + a_2,i(x-x_i)² + a_3,i(x-x_i)³ (1.1)

Относящегося к одному i-му звену сплайна, для которого x є [x_i,x_i+1], выбираются два дифференциальных условия согласования  с порядком Р₁={0; 2}, каждое из которых относится к концам отрезка [х_i, x_i+1]:

(1.2)

 

 

( 1.3)

 

 

Выбор значения p = 0 обусловлен тем, что отыскивается S₃(x) — интерполяционный многочлен, а выбор производных порядка р = 2 зависит от выбора способа. Условия (1.2), (1.3) задают параметры сплайна f_i,

M_i=f^’’(x_i), i= первые из которых (функциональные) являются определенными, а вторые (дифференциальные) - неопределенными.

Подстановка в соотношения (1.2), (1.3) полинома (1.1) приводит к системе четырех линейных алгебраических уравнений относительно неизвестных коэффициентов а_к,i (k = 0,1,2,3). Разрешая их, подставляя получен ные формулы для a_{k ,i} в (1.1) и распространяя этот результат на все частичные отрезки, имеем:

f

(1.4)

 

 

Где h_i+1 = x_i+1- x_i , Δf_i = f_i+1 – f_i, Δm_i = m_i+1 – m_i.

Таким образом, найдена общая формула  одного звена сплайна, которая выражается через определенные (f_i) и неопределенные (m_i) параметры

(i= ). Эти параметры обеспечивают непрерывность сплайна: 

, (i= ), и его вторых производных

, (i= ),  всех внутренних узлов.

Второй этап. При вычислении неопределенных параметров m_i(i= ) входящих в (1.4), с учетом условия стыковки для всех внутренних узлов устанавливается связь между f_i и mi. Эта связь получается путем записи правой части (1.4) для звена S_3,i-1(x) при хє[x_i-1,x_i], ее дифференцирования и приравнивания производной правой части (1.4) для звена S_3,i(x) при хє[x_i,x_i+1], соответствии с условием стыковки (р₂ = 1). Проведя несложные выкладки и распространив найденное соотношение на все внутренние узлы, получим параметрическую систему — трехдиагональную систему линейных алгебраических уравнений, устанавливающую линейную связь комбинаций определенных f_i и неопределенных параметров m_i, во внутренних узлах:

(1.5)

Данная система относительно m_i (i= ) является незамкнутой, так как не хватает двух уравнений. Для ее замыкания используют различные пути выбора аппроксимаций производных на концах (граничных условий):

1. В простейшем случае вторые производные на концах принимаются нулевыми, т.е.  m₀ = 0 ; m_n = 0.

 Такие условия называются  условиями натурального сплайна. 

2. Для первых двух и последних двух отрезков можно применить условие равенства третьей производной. Тогда получается соотношение

Δm_k-1/h_k=Δm_k/h_k+1 (k = 1,n- 1).

 

Отсюда следуют два трехточечных уравнения:

(1.6)

 

 

 

Соотношения (1.6) позволяют замкнуть систему (1.5) при h = var (неравномерная сетка). При h = const вторые производные на концах вычисляются по явным формулам второго порядка аппроксимации [6]:

(1.7)

 

 

Отметим, что как (1.6), так и формулы (1.7) для т₀ и т_п имеют порядки аппроксимации, соответствующие порядку сходимости S₃(x)  к f(x).

Третий этап. Определяются значения m_i (i= )  путем решения систем (1.5), (1.6) (параметрически прямая задача). Для этого используется метод прогонки, причем необходимо предварительно из первых двух и последних двух уравнений исключить слагаемые с т₂ и т_п-₂ соответственно. В случае равномерной сетки, как отмечено выше, система (1.5) замыкается аппроксимацией (1.7).

Четвертый этап. Формируются массивы коэффициентов a_0,i, a_1,i, a_2,i, a_3,i

для всех звеньев сплайна (i= ) и определяется глобальный интерполяционный сплайн:

путем составления многозвенной функции.

 

Пятый этап. При необходимости полученная сплайн-функция S₃(x) используется для вычисления значений функции и производных порядка р :

(p = 0,1,2…) в произвольных точках хє[x_i-1,x_i], (i= ) , или определенных интегралов:

Методика построения интерполяционного  кубического сплайна 

Используя систему (1.5) , условия т₀ = 0, т_п = О или (1.6) (при h = const можно использовать условия (1.7)), сформировать замкнутую систему относительно неопределенных параметров.

Вычислить коэффициенты системы (1.5), зависящие от шагов сетки Q„, и решить ее методом прогонки. В результате найти неопределенные параметры

m₀,m₁,…,m_n.

Записать выражения для звеньев  сплайна по формуле (1.4):

Вычислить коэффициенты звеньев и сформировать многозвенную сплайн-функцию

являющуюся сплайном дефекта q = 1.

 

 

 

 

3. Обзор методов решения СЛАУ с трехдиагональными матрицами методом фронтальной прогонки.

Метод применим в случае, когда матрица А — трехдиагональная. Общая постановка задачи имеет следующий вид.

Дана  система линейных алгебраических уравнений  с трехдиагональной матрицей А . Развернутая запись этой системы имеет вид

(1.8)

которому соответствует расширенная матрица

Здесь первое и последнее уравнения, содержащие по два слагаемых, могут рассматриваться  как краевые условия. Знак «—»  при коэффициенте β_i, взят для более удобного представления расчетных формул метода.

Требуется найти решение х_* =  (x_*1,…,x_*n)^T системы (1.8) методом исключения Гаусса.

Если  к (1.8) применить алгоритм прямого хода метода Гаусса, то вместо А₁ получится А₁:

Учитывая, что последний столбец в этой матрице соответствует правой части, и переходя к системе, включающей неизвестные, получаем рекуррентную формулу:

(1.9)

Соотношение (1.9) есть формула для обратного хода, а формулы для коэффициентов P_i,Q_i , которые называются прогоночными, определяются из (1.8), (1.9). Запишем (1.9) для индекса (i-1):

И подставим  в (1.4). Получим 

Приводя эту формулу к виду (1.9) и сравнивая полученное выражение с (1.9), получаем рекуррентные соотношения для P_i,Q_i :

(2.0)

 

Определение прогоночных коэффициентов по формулам (2.0) соответствует прямому ходу метода прогонки.

Обратный ход метода прогонки начинается с вычисления х_п. Для этого используется последнее уравнение, коэффициенты которого определены в прямом ходе, и последнее уравнение исходной системы:

Тогда определяется х_n:

Остальные значения неизвестных находятся рекуррентно по формуле (1.9).

Все соотношения  для выполнения вычислений получены. Тогда можно провести расчеты  по методу Гаусса, используя прямой и обратный ход.

 

 

Методика  решения задачи

Прямой ход.

1. Вычислить , (в (2.0) подставить a₁ =0).

     2.  Вычислить прогоночные коэффициенты P₂,Q₂; P₃,Q₃;…; P_n-1,Q_n-1; по формулам (2.0).

Обратный ход.

1. Найти

Значения  x_n-1, x_n-2,…, x₁, определить по формуле (1.9) :

Замечания.

1. Данный метод называется методом скалярной прогонки, так как при решении задачи на каждом i-ом шаге определяется скалярная величина  x_i (i= ).
2. Аналогичный подход используется для решения систем линейных алгебраических уравнений с пятидиагональными матрицами.
3. Алгоритм метода прогонки называется корректным, если для все

(i= )  , и если устойчивым, если |Р_i|<1, (i= ).

4. Достаточным условием коррекности и устойчивости прогонки является условие преобладания диагональных элементов в матрице А, в матрице a_i ≠ 0 и γ_i ≠ (i= ):

 

 

 

 

 

 

 

4.   Численный пример.

1. Равномерная сетка.

Пусть дана равномерная сетка, со значениями функции в точках.

Хi	F(xi)	ΔF(xi)	h
0	0	0.10017	0.05
0.05	0.10017	0.10117
0.1	0.20134	0.10318
0.15	0.30452

 

Мы так  же сразу можем найти ΔF(x_i), так как в дальнейших вычислениях эти данные нам пригодятся. ΔF(x_i) находится по формуле:

 

И найдем шаг сетки h. Шаг находится по формуле:

 

Запишем систему для внутренних узлов:

 

 

Так как  h = const, для нахождения значения выбираются формулы:

 

 

Подставляем в эти формулы данные с таблицы, и получаем следующие значения , .

 

 

 

 

 

Получается  система:

 

 

 

 

Далее эту  систему мы будем решать методом  прогонки (см. п. 1.2).

По методу прогонки нам надо вычислить:

; ; по формулам :

;  , где .

Найденные значения P: P₁ = 0; Р₂ = -0,25; Р₃ = -0,2667.

Найденные значения Q: Q₁ = 1.00169; Q₂ = 0.349578; Q₃= 1.193179; Q₄= 1.208.

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

 

.

Получаем  значения:  M₁ = 1.00169; M₂ = 0.131816; M₃ = 0.871046; M₄ = 1.208.

Все значения нам известны. Можно приступить к нахождения сплайна.

 

 

Для нахождения сплайна нам понадобится формула:

Где  h_i+1 =  x_i+1- x_i,  Δf_i  = f_i+1 - f_i , Δm_i = m_i+1 – m_i , i=0,1,2,3.

 

 

 Нам  нужно найти значение сплайна  в двух точках Х₁ = 0,06 и Х₂ = 0,11. Чтобы знать в какой из сплайнов подставлять Х, надо подставить в промежутки. Х₁ подходит для второго уравнения, Х₂ подходит для третьего. Подставляя Х мы находим значение сплайна в нужном Х.

Подставив 0,06 во второй сплайн, а 0,11в третий сплайн, нашли значение сплайнов, и они  равны:

S_3,1(0.06) = 0.120318, S_3,2 (0.11) = 0.221774

 

 

 

2. Неравномерная сетка.

Пусть дана неравномерная сетка, со значениями функции в точках Х.

I	Xi	F(Xi)	Δ F(Xi)	Hi
0	0.1	1.1052	0.0566	0.05
1	0.15	1.1618	0.0474	0.04
2	0.19	1.2092	0.0748	0.06
3	0.25	1.284

 

Δ F(X_i) = F(X_i+1) - F(X_i), h_i = h_i+1 – h_i , i=0,1,2,3.

Запишем систему для внутренних узлов:

 

 

Так как  h = var , для нахождения значения выбираются формулы:

 

 

 

Получается  система из 4х уравнений, подставим известные нам значения:

 

 

 

 

 

 

Эту систему  можно решить методом Гаусса:

20	-45	25	0	0
0,00833333	0,03	0,006667	0	0,053
0	0,006667	0,033333	0,01	0,061667
0	25	-8,33333	16,66667	0

 

1	-2,25	1,25	0	0
0	0,04875	-0,00375	0	0,053
0	0,006667	0,033333	0,01	0,061667
0	25	-8,33333	16,66667	0

 

1	-2,25	1,25	0	0
0	1	-0,07692	0	1,087179
0	0	0,033846	0,01	0,054419
0	0	-6,41026	16,66667	-27,1795

 

1	-2,25	1,25	0	0
0	1	-0,07692	0	1,087179
0	0	1	0,295455	1,607828
0	0	0	18,56061	-16,8729

 

1	-2,25	1,25	0	0
0	1	-0,07692	0	1,087179
0	0	1	0,295455	1,607828
0	0	0	1	-0,90907

 

По обратному  ходу найдем искомые M_i :

M₁ = 0.425397; M₂ = 1.231519; M₃ = 1.876417; M₄ = -0.90907;

Найдем ΔM_i = M_i+1 – M_i . M₁ = 0.806122; M₂ = 0.644898; M₃ = -2.78549.

Все значения нам известны. Можно приступить к нахождения сплайна.

Для нахождения сплайна нам понадобится формула:

Где  h_i+1 =  x_i+1- x_i,  Δf_i  = f_i+1 - f_i , Δm_i = m_i+1 – m_i , i=0,1,2,3.

 

 

 

 

 

Нам нужно  найти значение сплайна в двух точках Х₁ = 0,2 и Х₂ = 0,17. Чтобы знать в какой из сплайнов подставлять Х, надо подставить в промежутки. Х₁ подходит для второго уравнения, Х₂ подходит для третьего. Подставляя Х мы находим значение сплайна в нужном Х.

Подставив 0,17 во второй сплайн, а 0,2в третий сплайн, нашли значение сплайнов, и они  равны:

S_3,1(0.17) = 1.185168, S_3,2 (0.2) = 1.221476.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Глава II. Практическая реализация интерполирования кубическим сплайном

1. Реализация в С++.

Программа реализована с помощью трёх основных функций – В главной функции void main() реализована основная часть программы. Здесь хранится условие, по какому из методов будет решаться сплайн:

void main()

{

 clrscr();

 int ch;

 cout<<"Viberi setky:"<<endl;

 cout<<"1. Ravnomernaja"<<endl;

 cout<<"2.Neravnomernaja"<<endl;

 cin>>ch;

 if(ch==1) ravn();

 if(ch==2)neravn();

 cout<<"\n Vihod!";

   getch();

}

 

Функция void neravn() и void ravn(). Предназначены для вычисления всех необходимых данных для нахождения сплайна. К примеру часть программы которая вычисляет сплайн выглядит так :

cout<<"\n Vvedite zna4enie X - "<<endl;

   cin>>ot;

   if(ot>=x[0] && ot<x[1])

{

s3[0]=((1.0/h)*dfx[0]-(h/2)*ix[1]-(h/6)*dm[1])*(otx[0])+(ix[1]/2)*(pow(ot--x[0],2)) ;

s3[0]=s3[0]+(1.0/(6.0*h)*dm[1]*pow(ot-x[0],3)   ) ;

 s3[0]+=fx[0];

 cout<<"\n S3[0] = "<<s3[0];     }

 

  if(ot>=x[1] && ot<x[2])

{

s3[1]=((1.0/h)*dfx[1]-(h/2)*ix[2]-(h/6)*dm[2])*(ot-x[1])+(ix[2]/2)*(pow(ot--x[0],2)) ;

s3[1]=s3[1]+(1.0/(6.0*h)*dm[2]*pow(ot-x[1],3)   ) ;

s3[1]+=fx[1];

cout<<"\n S3[1] = "<<s3[1];  }

 

  if(ot>=x[2] && ot<x[3])

{

s3[2]=((1.0/h)*dfx[2]-(h/2)*ix[3]-(h/6)*dm[3])*(ot-x[2])+(ix[3]/2)*(pow(ot-x[1],2)) ;

s3[2]=s3[2]+(1.0/(6.0*h)*dm[2]*pow(ot-x[2],3)   ) ;

 s3[2]+=fx[2];

 cout<<"\n S3[2] = "<<s3[2]; }

 

2. Реализация в Exel.

Для вычисления интерполяционного сплайна требуется знать значения X и Y, которые будут записаны в сетку. В данном пакете реализованы вычисления как по равномерной сетке, так и по неравномерной сетке.

1.Равномерная сетка.

 

 

 

 

 

 

 

2.Неравномерная сетка.

 

 

 

 

 

3. Реализация в MathCad.

Для вычисления интерполяционного сплайна в  редакторе MathCad существует спефиальная функция.

-vx и vy, содержащие координаты x и y, через которые нужно провести кубический сплайн.

- vs := cspline( vx, vy). Вектор vs содержит вторые производные интерполяционной кривой в рассматриваемых точках. функция  cspline генерирует кривую сплайна, которая может быть кубическим полиномом в граничных точках.

- Чтобы найти  интерполируемое значение в произвольной  точке, скажем x0, вычислите interp(vs, vx,vy, x0), где vs, vx и vy — векторы, описанные ранее.

 

 

 

 

 

 

 

 

 

 

2.4 Сравнение полученных результатов.

Решая  данный интерполяционный сплайн   в Excel, C++, MathCad можно сделать вывод, что наиболее точными из трех вышеприведенных программ является  Excel. Погрешность незначительная, но всё же она есть (5-го знак после запятой).

Результаты вычислений в С++ по равномерной сетке:

 

Результаты вычислений в С++ по неравномерной сетке:

 

 

 

 

 

 

 

Результаты вычислений в MathCad:

Неравномерная сетка.

Равномерная сетка:

 

 

Результаты вычислений в Exсel:

 

Равномерная сетка:

 

  Неравномерная сетка:

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2.5 Блок-схема программной  реализации сплайн интерполяции.

 

Начало

Ввод: X_i,F(X_i)