Интерполяция функций формулами Стирлинга

«Интерполяция функций формулами Стирлинга»

Яковлева Н.С.

Г.Магнитогорск,2012

Содержание

Введение.............................................................................................................. 3

1. Задача.............................................................................................................. 4

2. Интерполяция функций................................................................................. 5

2.1 Интерполяционная формула Ньютона....................................................................................... 6

2.2 Интерполяционные формулы Гаусса......................................................................................... 7

2.3 Интерполяционная формула Стирлинга.................................................................................. 10

3. Разработка программного обеспечения..................................................... 11

4. Текст программы......................................................................................... 13

5. Тестирование программы............................................................................ 17

6. Результаты решения задачи........................................................................ 21

7. Построение графика функции на участке интерполирования................ 24

8. Заключение................................................................................................... 25

Литература........................................................................................................ 26

ВВЕДЕНИЕ

Целью моей курсовой работы является интерполяция функции, заданной таблично, методом Стирлинга (Бесселя) с использованием программного обеспечения, разработанного на любом языке программирования. Для этого я рассмотрела интерполяционные формулы Ньютона, Гаусса и, конечно же, Стирлинга. При их рассмотрении получаем, что полусумма двух интерполяционных формул Гаусса и даёт нам формулу Стирлинга. При разработке программное обеспечение, с целью автоматизации расчетов при выполнении интерполяции по формуле Стирлинга, разработана программа на языке Паскаль, в среде программирования Turbo Pascal 7.0. После того, как была установлена корректность выдаваемых программой результатов расчетов, я приступила к решению поставленной задачи. Представила протокол работы программы при вычислении значений функции на интервале интерполирования. Итогом моей курсовой стал график заданной мною функции. По данным, полученным с помощью программы, построили средствами табличного процесора Excel график функции на участке интерполирования.
1. ЗАДАЧА

Интерполировать функцию, заданную таблично, методом Стирлинга (Бесселя) с использованием программного обеспечения, разработанного на любом языке программирования.

Исходные данные:

i	xi	f(xi)
1	20	1002,3
2	40	541,7
3	60	116,87

2. Интерполяция функций

Под приближением функции f(x), заданной на интервале [a, b], будем понимать замену f(x) некоторой другой функцией P(x), близкой к исходной функции f(x). Простейшая задача, приводящая к приближению функций, состоит в следующем. Пусть на сетке {x_i} задана табличная функция f(x_i) = y_i. Требуется найти значения функции f(x) в точках x_j, не совпадающих с узлами исходной сетки x_i.

Часто приближающую функцию P(x) задают в следующей форме:

,

где  — система функций, заданных на [a, b] и являющихся гладкими (непрерывно дифференцируемыми); a_i — коэффициенты, которые выбирают таким образом, чтобы отклонение f(x) от P(x) было минимальным на заданном множестве X = {x}.

Если параметры a_i, i = 0, 1, …, n, определяются из условия совпадения значений f(x) и P(x) в узлах сетки x_i

f(x_i) = P(x_i),

то такой способ приближения называют интерполяцией, или интерполированием, а сетку {x_i} называют интерполяционной сеткой. При этом предполагается, что значения сетки {x_j}, в которой мы хотим вычислить P(x_j), не выходят за пределы границ сетки {x_i}:

a ≤ x_j ≤ b, j = 1, 2, …, m.

 

2.1 Интерполяционная формула Ньютона

Введем следующие понятия:

1.     Отношения  (i = 0, 1, …) называются разделенными разностями первого порядка для табличной функции f(x).

2.     Отношения  называются разделенными разностями второго порядка.

3.     Для разделенных разностей n-го порядка имеем

.

Лемма. Если f(x) = P(x) есть полином n-й степени, то его разделенная разность (n+1)-го порядка равна нулю, т. е.

для любой системы различных между собой чисел .

Пусть P(x) — полином степени n, такой, что P(x_i) = y_i = f(x_i) (i = 0, 1, …, n), f(x) — данная функция. Обозначим через P(x, x₀); P(x, x₀, x₁); …, P(x, x₀, x₁, …, x_n) разделенные разности полинома P(x), причем P(x, x₀, x₁, …, x_n) = 0.

По определению имеем

,

отсюда

.                                                                               (1)

Далее,

Подставляя это в (1), получим

Здесь последнее слагаемое равно нулю. Учитывая, что  — разделенные разности m-го порядка, которые совпадают с разделенными разностями для функции y = f(x) (т. к. P(x_i) = f(x_i)), то окончательно получим:

                                 (2)

 

2.2 Интерполяционные формулы Гаусса

Воспользуемся интерполяционной формулой Ньютона для неравных промежутков (2) и возьмем в качестве узлов  точки  Тогда

                (3)

Используя симметричность разделенных разностей относительно своих аргументов и их связь с конечными разностями , получим:

,                                                        (4)

.                                            (5)

Отсюда

         (6)

Обозначив

,

получим:

                   (7)

Это — интерполяционная формула Гаусса. В ней используются следующие разности (подчеркнуты черточкой):

x	f	f1	f2	f3	f4
x–3	f–3
x–2	f–2
x–1	f–1
x0
x1	f1
x2	f2
x3	f3

Если бы мы взяли узлы интерполирования в другом порядке, а именно: , то совершенно аналогично получили бы вторую формулу Гаусса:

                   (8)

Для того чтобы их можно было различить, будем называть первую из них интерполяционной формулой Гаусса для интерполирования вперед, а вторую — для интерполирования назад. Интерполяционная формула Гаусса для интерполирования назад использует следующие разности:

x	f	f1	f2	f3	f4
x–3	f–3
x–2	f–2
x–1	f–1
x0
x1	f1
x2	f2
x3	f3

 

2.3 Интерполяционная формула Стирлинга

Полусумма двух интерполяционных формул Гаусса даст нам:

                           (9)

так как

а

.

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

x	f	f1	f2	f3	f4
x–2	f–2
x–1	f–1
x0
x1	f1
x2	f2

Остаточный член формулы Стирлинга

Если последняя из использованных в формуле Стирлинга разностей имеет нечетный порядок, то остаточный член будет иметь вид

3.     Разработка программного обеспечения

С целью автоматизации расчетов при выполнении интерполяции по формуле Стирлинга разработана программа на языке Паскаль, в среде программирования Turbo Pascal 7.0.

Программа после запуска запрашивает у пользователя ввод:

·        координаты центрального узла (x₀);

·        шага таблицы (h);

·        количества точек в исходной таблице до и после центрального узла (n; общее число точек в таблице, таким образом, равно 2·n + 1).

После ввода этих данных программа рассчитывает координаты абсцисс точек таблицы и запрашивает у пользователя значения функции в этих точках. Когда исходные данных полностью определены, программа вычисляет коэффициенты  из (9) для i от 1 до 2·n; эти коэффициенты позволяют выполнять вычисления по формуле Стирлинга при произвольном значении аргумента. Далее программа предлагает пользователю вводить необходимые значения x, а для выхода из программы ввести в качестве x значение координаты центрального узла.

Блок-схема программы показана на рис. 1.

Рис. 1. Схема алгоритма (начало)

Рис. 1. Схема алгоритма (продолжение)

Рис. 1. Схема алгоритма (окончание)

4.     Текст программы

program Stirling;{интерполяция табличных данных по формуле Стирлинга}

const NN = 50;{предельное число отсчетов от центрального узла вперед и назад (общее число узлов 2 * NN + 1)}

type ar = array [-NN..NN] of double;

var x, f, {исходные табличные данные}

    df: ar; {конечные разности}

    pf: array [1..2*NN] of double; {коэффициенты интерполяционного многочлена}

    i, k, n, dn, j: integer;

    h, t, p, fc, t2: double;

    ct: array [0..1] of double;

BEGIN

     writeln('Интерполяция по формуле Стирлинга');

     write('Введите координату центрального узла (x0): ');

     readln(x[0]);

     write('Введите шаг таблицы (h): ');

     readln(h);

     write('Введите кол-во узлов по каждую сторону ', 'от центрального (n): ');

     readln(n);

     dn := 2 * n; {сохраняем 2*n, чтобы не вычислять

        это значение каждый раз}

     writeln('Введите исходные данные:');

     for i := -n to n do

         begin

         x[i] := x[0] + i * h;

         write('f(', x[i]:0:2, ') = ');

         readln(f[i]);

         df[i] := f[i]; {заполняем массив конечных разностей значениями функции}

         end;

     for i := 1 to dn d {расчет коэффициентов

           интерполяционного многочлена}

         begin

         for k := -n to n - i do

             df[k] := df[k + 1] - df[k]; {вычисление

                   конечных разностей i-го порядка}

         if odd(i) then

            begin

            j := (dn + 1 - i) div 2 - n;

            pf[i] := (df[j - 1] + df[j]) / 2; {полусумма}

            end

         else pf[i] := df[(dn - i) div 2 - n];

         end;

     writeln('Вводите x для интерполирования');

     writeln('(для выхода из программы введите x = ',

                   x[0]:0:2, ')');

     repeat

           write('x = ');

           readln(t);

           t := (t - x[0]) / h;

           t2 := sqr(t); {сохраняем, чтобы не вычислять

              каждый раз}

           p := f[0]; {первый член интерполяционной формулы}

           fc := 1; {для вычисления факториала}

           ct[1] := t; {начальное значение аргумента для нечетных членов интерполяционной формулы}

           ct[0] := t2; {начальное значение аргумента для четных членов интерполяционной формулы}

           for i := 1 to dn do {вычисления по формуле}

               begin

               fc := fc * i; {факториал}

               k := i mod 2; {k = 0 для четных членов, и 1 - для нечетных}

               if i > 2 then{умножать начальные значения аргументов на (t^2 - m^2) нужно только для членов после второго}

                  ct[k] := ct[k] * (t2 - sqr((i - 1) div 2));

               p := p + ct[k] * pf[i] / fc;

               end;

           writeln('f(', x[0] + t * h:0:2, ') = ', p:0:8);

     until t = 0;

END.

5.     Тестирование программы

В [3, с. 128] приведены результаты расчета значения sin(14º) по известным значениям синуса при углах 9º, 12º, 15º, 18º, 21º. Выполним тот же расчет с помощью разработанной программы:

Интерполяция по формуле Стирлинга

Введите координату центрального узла (x0): 15

Введите шаг таблицы (h): 3

Введите кол-во узлов по каждую сторону от центрального (n): 2

Введите исходные данные:

f(9.00) = 0.156434

f(12.00) = 0.207912

f(15.00) = 0.258819

f(18.00) = 0.309017

f(21.00) = 0.358368

Вводите x для интерполирования

(для выхода из программы введите x = 15.00)

x = 14

f(14.00) = 0.24192196

x = 15

f(15.00) = 0.25881900

Полученные результаты согласуются с результатами, приведенными в [1], с точностью до 6 знаков после запятой — именно с такой точностью расчет проделан в первоисточнике, причем указано, что все найденные знаки в результате — верные.

6.     Результаты решения задачи

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

Интерполяция по формуле Стирлинга

Введите координату центрального узла (x0): 40

Введите шаг таблицы (h): 20

Введите кол-во узлов по каждую сторону от центрального (n): 1

Введите исходные данные:

f(20.00) = 1002.3

f(40.00) = 541.7

f(60.00) = 116.87

Вводите x для интерполирования

(для выхода из программы введите x = 40.00)

x = 20

f(20.00) = 1002.30000000

x = 22

f(22.00) = 954.63035000

x = 24

f(24.00) = 907.31840000

x = 26

f(26.00) = 860.36415000

x = 28

f(28.00) = 813.76760000

x = 30

f(30.00) = 767.52875000

x = 32

f(32.00) = 721.64760000

x = 34

f(34.00) = 676.12415000

x = 36

f(36.00) = 630.95840000

x = 38

f(38.00) = 586.15035000

x = 42

f(42.00) = 497.60735000

x = 44

f(44.00) = 453.87240000

x = 46

f(46.00) = 410.49515000

x = 48

f(48.00) = 367.47560000

x = 50

f(50.00) = 324.81375000

x = 52

f(52.00) = 282.50960000

x = 54

f(54.00) = 240.56315000

x = 56

f(56.00) = 198.97440000

x = 58

f(58.00) = 157.74335000

x = 60

f(60.00) = 116.87000000

x = 40

f(40.00) = 541.70000000

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

7.     Построение графика функции на участке интерполирования

По данным, полученным с помощью программы, построим средствами табличного процессора Excel график функции на участке интерполирования. График представлен на рис. 2.

Рис. 2. График функции на участке интерполирования

ЗАКЛЮЧЕНИЕ

В своей работе я рассмотрела интерполяционную формулу Стирлинга. Для полученной формулы составила программу на языке Паскаль, в среде программирования Turbo Pascal 7.0. Для проверки корректности выводимых ею результатов в [3, с. 128] взяла уже известные значения при углах 9º, 12º, 15º, 18º, 21º sin(14º). Выполнив те же расчеты с помощью разработанной программы, установила, что она работает исправно. После того, как была установлена корректность выдаваемых программой результатов, я приступила к решению поставленной задачи. Итогом моей курсовой стал график заданной мною функции. По данным, полученным с помощью программы, построили посредствами табличного процессора Excel график функции на участке интерполирования.

Литература

1.     Амосов А. А., Дубинский Ю. А., Копченова Н. В. Вычислительные методы для инженеров: Учебное пособие. — М.: Высшая школа, 1994. — 544 с.

2.     Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы. — М.: Наука, 2001. — 630 с.

3.     Березин И. С., Жидков Н. П. Методы вычислений. Т. 1. — М.: Гос. изд-во физ.-мат. литературы, 1962. — 464 с.

4.     Демидович Б. П., Марон И. А. Основы вычислительной математики. — М.: Издательство ФМЛ, 1960. — 659 с.

5.     Демидович Б. П., Марон И. А., Шувалова Э. З. Численные методы анализа. — М.: Издательство ФМЛ, 1963. — 400 с.

6.     Епанешников А. М., Епанешников В. А. Программирование в среде Turbo Pascal 7.0. — М.: “ДИАЛОГ-МИФИ”, 1993. — 288 с.

7.     Охлопков Н. М. Численные методы и вычислительные алгоритмы. Часть 1: Учебное пособие. — Якутск: Издательство ЯГУ, 1994. — 108 с.

8.     Фаронов В. В. Турбо Паскаль 7.0. Начальный курс. Учебное пособие — М.: Нолидж, 1997. — 616 с.