Интерполяционные кубические сплайны

Введение

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

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

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

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

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

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

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

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