Интерполяция функций формулами Стирлинга
«Интерполяция функций формулами Стирлинга»
Яковлева Н.С.
Г.Магнитогорск,2012
Содержание13 TOC \o "1-3" \u 14
Введение 3
1. Задача 4
2. Интерполяция функций 13 PAGEREF _Toc325633450 \h 14515
2.1 Интерполяционная формула Ньютона 6
2.2 Интерполяционные формулы Гаусса 7
2.3 Интерполяционная формула Стирлинга 10
3. Разработка программного обеспечения 13 PAGEREF _Toc325633454 \h 141115
4. Текст программы 13
5. Тестирование программы 13 PAGEREF _Toc325633456 \h 141715
6. Результаты решения задачи 21
7. Построение графика функции на участке интерполирования 24
8. Заключение 25
Литература 26
15 ВВЕДЕНИЕ
Целью моей курсовой работы является интерполяция функции, заданной таблично, методом Стирлинга (Бесселя) с использованием программного обеспечения, разработанного на любом языке программирования. Для этого я рассмотрела интерполяционные формулы Ньютона, Гаусса и, конечно же, Стирлинга. При их рассмотрении получаем, что полусумма двух интерполяционных формул Гаусса и даёт нам формулу Стирлинга. При разработке программное обеспечение, с целью автоматизации расчетов при выполнении интерполяции по формуле Стирлинга, разработана программа на языке Паскаль, в среде программирования 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). Простейшая задача, приводящая к приближению функций, состоит в следующем. Пусть на сетке {xi} задана табличная функция f(xi) = yi. Требуется найти значения функции f(x) в точках xj, не совпадающих с узлами исходной сетки xi.
Часто приближающую функцию P(x) задают в следующей форме:
13 EMBED Equation.3 1415,
где 13 EMBED Equation.3 1415 система функций, заданных на [a, b] и являющихся гладкими (непрерывно дифференцируемыми); ai коэффициенты, которые выбирают таким образом, чтобы отклонение f(x) от P(x) было минимальным на заданном множестве X = {x}.
Если параметры ai, i = 0, 1, , n, определяются из условия совпадения значений f(x) и P(x) в узлах сетки xi
f(xi) = P(xi),
то такой способ приближения называют интерполяцией, или интерполированием, а сетку {xi} называют интерполяционной сеткой. При этом предполагается, что значения сетки {xj}, в которой мы хотим вычислить P(xj), не выходят за пределы границ сетки {xi}:
a
· xj
· b, j = 1, 2, , m.
2.1 Интерполяционная формула Ньютона
Введем следующие понятия:
Отношения 13 EMBED Equation.3 1415 (i = 0, 1, ) называются разделенными разностями первого порядка для табличной функции f(x).
Отношения 13 EMBED Equation.3 1415 называются разделенными разностями второго порядка.
Для разделенных разностей n-го порядка имеем
13 EMBED Equation.3 1415.
Лемма. Если f(x) = P(x) есть полином n-й степени, то его разделенная разность (n+1)-го порядка равна нулю, т. е.
13 EMBED Equation.3 1415
для любой системы различных между собой чисел 13 EMBED Equation.3 1415.
Пусть P(x) полином степени n, такой, что P(xi) = yi = f(xi) (i = 0, 1, , n), f(x) данная функция. Обозначим через P(x, x0); P(x, x0, x1); , P(x, x0, x1, , xn) разделенные разности полинома P(x), причем P(x, x0, x1, , xn) = 0.
По определению имеем
13 EMBED Equation.3 1415,
отсюда
13 EMBED Equation.3 1415. (1)
Далее,
13 EMBED Equation.3 1415
13 EMBED Equation.3 1415
13 EMBED Equation.3 1415
Подставляя это в (1), получим
13 EMBED Equation.3 1415
Здесь последнее слагаемое равно нулю. Учитывая, что 13 EMBED Equation.3 1415 разделенные разности m-го порядка, которые совпадают с разделенными разностями для функции y = f(x) (т. к. P(xi) = f(xi)), то окончательно получим:
13 EMBED Equation.3 1415 (2)
2.2 Интерполяционные формулы Гаусса
Воспользуемся интерполяционной формулой Ньютона для неравных промежутков (2) и возьмем в качестве узлов 13 EMBED Equation.3 1415 точки 13 EMBED Equation.3 1415 Тогда
13 EMBED Equation.3 1415 (3)
Используя симметричность разделенных разностей относительно своих аргументов и их связь с конечными разностями 13 EMBED Equation.3 1415, получим:
13 EMBED Equation.3 1415, (4)
13 EMBED Equation.3 1415. (5)
Отсюда
13 EMBED Equation.3 1415 (6)
Обозначив
13 EMBED Equation.3 1415,
получим:
13 EMBED Equation.3 1415 (7)
Это интерполяционная формула Гаусса. В ней используются следующие разности (подчеркнуты черточкой):
x
f
f1
f2
f3
f4
x–3
f–3
13 EMBED Equation.3 1415
x–2
f–2
13 EMBED Equation.3 1415
13 EMBED Equation.3 1415
13 EMBED Equation.3 1415
x–1
f–1
13 EMBED Equation.3 1415
·–
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·–
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·–
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·–
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·–
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·–
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·–
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·Если бы мы взяли узлы интерполирования в другом порядке, а именно: 13 EMBED Equation.3 1415, то совершенно аналогично получили бы вторую формулу Гаусса:
13 EMBED Equation.3 1415 (8)
Для того чтобы их можно было различить, будем называть первую из них интерполяционной формулой Гаусса для интерполирования вперед, а вторую для интерполирования назад. Интерполяционная формула Гаусса для интерполирования назад использует следующие разности:
x
f
f1
f2
f3
f4
x–3
f–3
13 EMBED Equation.3 1415
x–2
f–2
13 E
·
·
·
·
·
·
·
·
·
·
·
·–
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·–
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·–
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·–
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·–
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·–
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·–
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·–
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·–
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·
·Интерполяционная формула Стирлинга
Полусумма двух интерполяционных формул Гаусса даст нам:
13 EMBED Equation.3 1415 (9)
так как
13 EMBED Equation.3 1415
а
13 EMBED Equation.3 1415.
Мы получили формулу Стирлинга. В ней используются разности четного порядка с индексом 0 и полусуммы разностей нечетного порядка с индексами 13 EMBED Equation.3 1415 и 13 EMBED Equation.3 1415, как это показано в следующей таблице:
x
f
f1
f2
f3
f4
x–2
f–2
13 EMBED Equation.3 1415
x–1
f–1
13 EMBED Equation.3 1415
13 EMBED Equation.3 1415
13 EMBED Equation.3 1415
x0
13 EMBED Equation.3 1415
13 EMBED Equation.3 1415
13 EMBED Equation.3 1415
x1
f1
13 EMBED Equation.3 1415
13 EMBED Equation.3 1415
x2
f2
Остаточный член формулы Стирлинга
13 EMBED Equation.3 1415
Если последняя из использованных в формуле Стирлинга разностей имеет нечетный порядок, то остаточный член будет иметь вид
13 EMBED Equation.3 1415
Разработка программного обеспечения
С целью автоматизации расчетов при выполнении интерполяции по формуле Стирлинга разработана программа на языке Паскаль, в среде программирования Turbo Pascal 7.0.
Программа после запуска запрашивает у пользователя ввод:
координаты центрального узла (x0);
шага таблицы (h);
количества точек в исходной таблице до и после центрального узла (n; общее число точек в таблице, таким образом, равно 2·n + 1).
После ввода этих данных программа рассчитывает координаты абсцисс точек таблицы и запрашивает у пользователя значения функции в этих точках. Когда исходные данных полностью определены, программа вычисляет коэффициенты 13 EMBED Equation.3 1415 из (9) для i от 1 до 2·n; эти коэффициенты позволяют выполнять вычисления по формуле Стирлинга при произвольном значении аргумента. Далее программа предлагает пользователю вводить необходимые значения x, а для выхода из программы ввести в качестве x значение координаты центрального узла.
Блок-схема программы показана на рис. 1.
Рис. 13 SEQ Рис. \* ARABIC 14115. Схема алгоритма (начало)
Рис. 1. Схема алгоритма (продолжение)
Рис. 1. Схема алгоритма (окончание)
Текст программы
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.
Тестирование программы
В [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 знаков после запятой именно с такой точностью расчет проделан в первоисточнике, причем указано, что все найденные знаки в результате верные.
Результаты решения задачи
После того, как была установлена корректность выдаваемых программой результатов расчетов, можно приступить к решению поставленной задачи. Ниже представлен протокол работы программы при вычислении значений функции на интервале интерполирования с шагом 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
Как показывают результаты расчетов, в узлах интерполяции программа выдает значения, в точности равные табличным это еще раз позволяет убедиться в правильности ее работы.
Построение графика функции на участке интерполирования
По данным, полученным с помощью программы, построим средствами табличного процессора Excel график функции на участке интерполирования. График представлен на рис. 2.
13 EMBED Excel.Chart.8 \s 1415
Рис. 13 SEQ Рис. \* ARABIC 14215. График функции на участке интерполирования
ЗАКЛЮЧЕНИЕ
В своей работе я рассмотрела интерполяционную формулу Стирлинга. Для полученной формулы составила программу на языке Паскаль, в среде программирования Turbo Pascal 7.0. Для проверки корректности выводимых ею результатов в [3, с. 128] взяла уже известные значения при углах 9є, 12є, 15є, 18є, 21є sin(14є). Выполнив те же расчеты с помощью разработанной программы, установила, что она работает исправно. После того, как была установлена корректность выдаваемых программой результатов, я приступила к решению поставленной задачи. Итогом моей курсовой стал график заданной мною функции. По данным, полученным с помощью программы, построили посредствами табличного процессора Excel график функции на участке интерполирования.
Литература
Амосов А. А., Дубинский Ю. А., Копченова Н. В. Вычислительные методы для инженеров: Учебное пособие. М.: Высшая школа, 1994. 544 с.
Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы. М.: Наука, 2001. 630 с.
Березин И. С., Жидков Н. П. Методы вычислений. Т. 1. М.: Гос. изд-во физ.-мат. литературы, 1962. 464 с.
Демидович Б. П., Марон И. А. Основы вычислительной математики. М.: Издательство ФМЛ, 1960. 659 с.
Демидович Б. П., Марон И. А., Шувалова Э. З. Численные методы анализа. М.: Издательство ФМЛ, 1963. 400 с.
Епанешников А. М., Епанешников В. А. Программирование в среде Turbo Pascal 7.0. М.: “ДИАЛОГ-МИФИ”, 1993. 288 с.
Охлопков Н. М. Численные методы и вычислительные алгоритмы. Часть 1: Учебное пособие. Якутск: Издательство ЯГУ, 1994. 108 с.
Фаронов В. В. Турбо Паскаль 7.0. Начальный курс. Учебное пособие М.: Нолидж, 1997. 616 с.
13PAGE 15
13PAGE \* MERGEFORMAT142615