Содержание
Введение4
1 Понятие интерполяции5
2 Интерполяционные формулы Ньютона6
2.1 Первая интерполяционная формула Ньютона9
2.2 Вторая интерполяционная формула Ньютона14
3 Формула Лагранжа20
3.1 Описание пользовательского интерфейса22
4 Формула Эткена25
5 Интерполяционная формула Стирлинга26
6 Интерполяционная формула Бесселя:27
Литература28
Приложение 129
Приложение 240
Приложение 3.46
Выдержка из текста работы
|
|
Введение |
2 |
|
1. |
Интерполирование функций полиномом Лагранжа. |
|
|
1.1. |
Теоретические основы метода |
3 |
|
1.2. |
Синтез полинома Лагранжа для заданной функции. Оценка точности интерполяции. |
4 |
|
|
Выводы |
11 |
|
2. |
Аппроксимация функций полиномом Ньютона. |
|
|
2.1. |
Теоретические основы метода. |
11 |
|
2.2. |
Синтез полинома Ньютона |
12 |
|
2.2.1. |
Интерполяция «вперед» |
12 |
|
2.2.2. |
Интерполяция «назад» |
19 |
|
|
Выводы |
19 |
|
3. |
Интерполяция функций методом наименьших квадратов.
|
|
|
3.1. |
Теоретические основы метода |
20 |
|
3.2. |
Составление аппроксимирующего полинома для выбранной функции и оценка точности аппроксимации |
21 |
|
3.2.1. |
Линейная аппроксимация |
22 |
|
3.2.2. |
Квадратичная аппроксимация |
24 |
|
|
Выводы |
34 |
|
|
Заключение |
35 |
|
|
Список литературы |
36 |
Введение
В задачах теории колебаний, электродинамики и других разделах прикладной математики широко применяется аппроксимация функций при описании физических параметров сред. В задачах вычислительной математики аппроксимация функций является основой для разработки многих методов и алгоритмов.
Аппроксимация, или приближение — научный метод, состоящий в замене одних объектов другими, в том или ином смысле близкими к исходным, но более простыми. Иными словами, аппроксимация некоторой функции y=f(x) заключается в замене другой функцией g(x, a0,a1,…, an) так, чтобы отклонение g от f(x) удовлетворяло в некоторой области (на множестве Х) определенному условию.
При этом функция g(x, a0,a1,…, an) обычно выбирается с учетом специфических особенностей рассматриваемой функции f(x). Если множество Х дискретно, то приближение называется точечным, если же Х есть отрезок a?x?b, то приближение называется интегральным.
Частным случаем аппроксимации является интерполяция.
Обычно задача аппроксимации распадается на две части:
1. Сначала устанавливают вид зависимости y=f(x) и, соответственно вид эмпирической формулы, то есть решают, является ли она линейной, квадратичной, логарифмической или какой-либо другой.
2. После этого определяются численные значения неизвестных параметров выбранной эмпирической формулы, для которых приближение к заданной функции оказывается наилучшим.
Иной подход к решению задачи аппроксимации:
— Если нет каких-либо теоретических соображений для подбора вида формулы, обычно выбирают функциональную зависимость из числа наиболее простых, сравнивая их графики с графиком заданной функции.
— После выбора вида формулы определяют ее параметры. Для наилучшего выбора параметров задают меру близости аппроксимации экспериментальных данных. Во многих случаях, в особенности если функция f(x) задана графиком или таблицей (на дискретном множестве точек), для оценки степени приближения рассматривают разности f(x)?(xi) для точек x0, x1, …, xn.
Для практики важен случай аппроксимации функции многочленами, т.е.
F(x)=a0+a1x+a2x2+…+anxn (1)
1 Интерполяционный полином Лагранжа.
1.1. Теоретические основы метода.
Интерполяционный полином Лагранжа — многочлен минимальной степени, принимающий данные значения в данном наборе точек. Для n+1 пар чисел (x0,y0),(x1,y1)… (xn,yn), где все xj различны, существует единственный многочлен L(x) степени не более n, для которого L(xj)=yj.
Интерполяционный полином Лагранжа имеет вид:
где li(x) – базисные полиномы, определяющиеся по формуле:
li(x) обладают следующими свойствами:
— являются многочленами степени n;
— li(xi) =1;
— li(xj) =0 при j?i.
Отсюда следует, что L(x), как линейная комбинация li(x), может иметь степень не больше n, и L(xi)=yi.
Оценку погрешности определяют по формуле:
|f(x)-Ln(x)|?|Ln+1(x)-Ln(x)| (4)
Задание №1:
Выполнить интерполяцию таблично заданной функции y=f(x).
|
i |
Х |
Y |
|
0 |
0.0 |
0.979498 |
|
1 |
0.1 |
1.060831 |
|
2 |
0.2 |
1.138837 |
|
3 |
0.3 |
1.213312 |
|
4 |
0.4 |
1.284076 |
|
5 |
0.5 |
1.350965 |
|
6 |
0.6 |
1.413842 |
|
7 |
0.7 |
1.472590 |
|
8 |
0.8 |
1.527116 |
|
9 |
0.9 |
1.577351 |
Вычислить значение интерполирующей функции в точке a=0.865. Оценить погрешность интерполяции.
1.2 Синтез полинома Лагранжа.
Перенумеруем узлы следующим образом: x0 – ближайший к заданной точке a узел и т.д.
|
i |
Х |
Y |
|
0 |
0.9 |
1.577351 |
|
1 |
0.8 |
1.527116 |
|
2 |
0.7 |
1.472590 |
|
3 |
0.6 |
1.413842 |
|
4 |
0.5 |
1.350965 |
|
5 |
0.4 |
1.284076 |
|
6 |
0.3 |
1.213312 |
|
7 |
0.2 |
1.138837 |
|
8 |
0.1 |
1.060831 |
|
9 |
0.0 |
0.979498 |
Построим интерполяционный многочлен Лагранжа для n=2. В этом случае он будет иметь вид:
L2(x)=l0*y0+l1*y1+l2*y2.
l0= ;
l1= ;
l2= ;
При увеличении порядка многочлена до n=4, полином Лагранжа будет выглядеть следующим образом:
L2(x)=l0*y0+l1*y1+l2*y2+l3*y3+l4*y4.
При этом, все базисные полиномы необходимо полностью пересчитать:
l0= ;
l1= ;
l2= ;
l3= ;
l4= ;
И, наконец, составим полином Лагранжа при n=9:
L2(x)=l0*y0+l1*y1+l2*y2+l3*y3+l4*y4+l5*y5+l6*y6+l7*y7+l8*y8+l9*y9.
Вновь пересчитываем базисные полиномы:
l0= ;
l1= ;
l2= ;
l3= ;
l4= ;
l5= ;
l6= ;
l7= ;
l8= ;
l9= .
Интерполяция полиномом Лагранжа возможна для произвольно заданных узлов интерполяции. Однако из приведенных выше расчетов видно, что добавление хотя бы одного нового узла интерполяции влечет за собой не только появление нового слагаемого, но и пересчет всех ранее вычисленных слагаемых.
Блок схема алгоритма интерполяции полиномом Лагранжа представлена на рис.1.
Входные данные:
— количество узлов интерполяции N;
— значение аргумента х в заданной точке;
— массив известных значений аргументов х(i) и функции y(i), сведенных в один двумерный массив arg(i,j). При этом, в связи с тем что значения аргументов изменяются в заданном диапазоне с равным интервалом, их значения вычисляются автоматически и записываются в первый столбец массива arg(i,1). Значения функции вводятся с клавиатуры и сохраняются во втором столбце (arg(i,2));
— порядок многочлена Лагранжа k.
Выходные данные:
— массив табличных значений аргументов arg(i,1) и функции arg(i,2);
— массив значений многочленов Лагранжа L();
— погрешность вычислений Е.
Программа включает в себя функцию lagrange блок-схема алгоритма которой, представлена на рис. 2. Функция принимает входные данные: порядок полинома Лагранжа, значение аргумента в заданной точке, массив табличных значений аргументов и функции в виде двумерного массива. И возвращает значение многочлена Лагранжа порядка k в точке х.
|
|
|
L(i)=lagrange |
|
i=0,N |
|
E=L(k)-L(k-1) |
|
flag=i |
|
Начало |
|
Ввод N, X |
|
Ввод x(i),y(i) |
|
x=x(i) |
|
flag=0 |
|
Ввод k |
|
i=1,k |
|
i=1,n |
|
Вывод x(i),y(i) |
|
Вывод x, L(k), E |
|
Вывод x,y(flag) |
|
Конец |
|
Рис. 1 Блок-схема алгоритма вычисления многочлена Лагранжа |
|
lagrange(k,x,arg) |
|
i=0,k |
|
l=arg(i,2) |
|
L=0 |
|
j=0,k |
|
i<>j |
|
l= |
|
L=L+l |
|
return L |
|
Рис.2 Блок схема алгоритма функции lagrange |
Текст программы вычисления многочлена Лагранжа.
// Method_L.cpp : main project file.
#include «stdafx.h»
#include <iostream>
#include<conio.h>
#include<cmath>
#include<iomanip>
#define max 20
using namespace System;
using namespace std;
double lagrange(int k,double x, double arg[][max]) /*объявление фнукции вычисления полинома Лагранжа, где «к»-колич-во узлов интерполяции, х-значение аргумента «х», arg- двумерный массив табличных значений аргументов и функции, в котором arg[i][1] – массив значений аргументов, arg[i][2] – массив значений функции*/
double l,L;
L=0;
for(int i=0;i<=k;i++)
l=arg[i][2];
for(int j=0;j<=k;j++)
if(i!=j)
l=(x-arg[j][1])/(arg[i][1]-arg[j][1])*l; //вычисление значения слагаемого полинома, при условии что «i»<>»j»
L=L+l; //подсчет суммы слагаемых полинома
return L; //возврат результата вычислений полинома вызывающей функции
int main() //начало основной программы
double arg[max][max],x,L[max],E;
float h;
int n,k,flag;
cout<<«Quantity of entered values N=»;
cin>>n; //определение количества табличных значений функции
cout<<«\n»;
cout<<«Enter value of argument X=»;
cin>>x; //ввод значения «х»
cout<<endl;
h=0.1;
flag=0;
for (int i=0;i<n;i++)
arg[i][1]=i*h; //заполнение массива табличных значений аргумента функции «x[] «
cout<<«Enter values of function Y»<<i<<«= «;
cin>>arg[i][2]; //заполнение массива табличных значений функции «y[]»
cout<<«\n»;
if(arg[i][1]==x) //проверка равенства введенного значения аргумента «х» со зна-
flag=i; чением из массива табличных аргументов «х[]»
cout<<«Enter quantity of knots of interpolation K=»;
cin>>k; //определение количества узлов интерполяции
cout<<«\n»;
if (flag==0) //если значение введенного аргумента «х» не равно табличному значению, то вычисляется полином Лагранжа
{
for (int i=1;i<=k;i++)
L[i]=lagrange(i,x,arg);
E=fabs(L[k]-L[k-1]); //вычисление погрешности
cout<<setw(5)<<«X(i)»<<setw(10)<<«Y(i)»<<endl;
for(int i=0;i<n;i++)
cout<<setw(5)<<arg[i][1]<<setw(10)<<arg[i][2]<<«\n»; //вывод на экран табличных значений аргумента и функции
cout<<endl;
cout<<«X=»<<x<<» L(«<<k<<«)=»<<L[k]<<» E=»<<E<<endl; //вывод «х», вычисленных значений функции полиномом Лагранжа и погрешности
}
else
{
cout<<«X=»<<x<<» Y=»<<arg[flag][2]<<endl; //вывод значения аргумента «х» и функции «y», в случае совпадения значения «х» с табличным
}
getch();
return 0;
Результаты интерполяции заданной функции полиномом Лагранжа различного порядка сведены в таблицу 1. На рисунке 3 представлен график исходной функции. На рисунке 4 представлен график функции, интерполированной полиномом Лагранжа различного порядка, в окрестности точки х=0.865.
Таблица 1.
|
Х |
n |
Ln(x) |
E |
|
0.865 |
2 |
1.572950 |
0.11 |
|
0.865 |
5 |
1.559700 |
7.7*10-4 |
|
0.865 |
9 |
1.560260 |
3.6*10-5 |
|
Рис. 3. График заданной функции. |
|
Х |
|
Y |
|
Рис. 4. График функции, интерполированной полиномом Лагранжа различного порядка, в окрестности точки х=0.865. |
|
Х |
|
Y |
Как видно из графика, изображенного на рис. 4, использование для интерполяции заданной функции полинома Лагранжа второго порядка, дает большую погрешность и явное отклонение графика функции в заданной точке от общего графика (рис.3). Тогда как при использовании полинома пятого порядка, точность вычислений значительно возрастает.
Выводы:
1. Формула Лагранжа может быть использована для произвольно заданных узлов интерполяции.
2. Добавление к уже взятым узлам интерполяции хотя бы одного нового, влечет за собой не только добавление нового базисного полинома, но и необходимость пересчета уже посчитанных базисных полиномов.
2. Аппроксимация функций полиномом Ньютона.
1.1. Теоретические основы метода.
Аппроксимация по формулам Ньютона используется если интерполируемая функция f(x) задана в (n+1) равноотстоящих узлах, т.е. xi+1—xi=?xi=h=const.
Общий вид полинома Ньютона:
Pn(x)=a0+a1(x-x0)+a2(x-x0)(x-x1)+…+an(x-x0)(x-x1)*…*(x-xn-1). (5)
Из условия интерполяции:
yi=f(xi)=Pn(xi), i= следует определение коэффициентов полинома Ньютона через конечные разности:
ai= , i=0,1,2,…n. (6)
В общем случае (для всех значений функции, определенных в узлах) конечные разности k-го порядка имеют вид:
?kyi= ?k-1yi+1— ?k-1yi (7)
Определим конечные разности:
?y0=y1-y0, ?y1=y2-y1, … ?yi=yi+1-yi,
?2y0=y1-y0, ?2y1=y2-y1, … ?2yi=yi+1-yi.
Найденные конечные разности записываются в виде таблицы:
|
x |
y |
?y |
?2y |
?3y |
|
|
x0 |
y0 |
?y0 |
?2y0 |
?3y0 |
|
|
x1 |
y1 |
?y1 |
?2y1 |
?3y1 |
|
|
x2 |
y2 |
?y2 |
?2y2 |
?3y2 |
|
|
|
|
|
|
|
|
Первая формула Ньютона имеет вид:
Pn(x)=y0+q ?y0+ ?2y0+…+ ?ny0 (8)
q= , где x0 – ближайший к точке х узел слева.
Если таблица значений функции конечна, то наивысший порядок интерполяционного многочлена равен n, где (n+1) – количество узлов интерполяции. Если длина таблицы задания функции не ограничена, то порядок интерполирования n определяется тем, что каждая конечная разность или их сумма не должна превышать некоторой постоянной величины Е:
Формула (8) применяется для интерполяции в начале таблицы («интерполяция вперед»), т.к. в этом случае имеется возможность, взяв большее число узлов, увеличить точность интерполирования.
Для интерполяции в конце таблицы («интерполяция назад») используется вторая интерполяционная формула Ньютона:
Pn(x)=yn+q ?yn-1+ ?2yn-2+…+ ?ny0 (11)
q= , где xn – ближайший к точке х узел справа.
Оценки погрешностей формул (8) и (11) имеют соответственно вид (12) и (13).
Rn(x)
Rn(x)
2.2. Синтез полинома Ньютона.
2.2.1.Интерполяция «вперед».
По формуле (7) найдем конечные разности и занесем их в таблицу 2.
Таблица 2.
|
Х |
Y |
?y |
?2y |
?3y |
?4y |
?5y |
?6y |
?7y |
?8y |
?9y |
|
0 |
0,979498 |
0,081333 |
-0,00333 |
-0,0002 |
2,4E-05 |
-8E-06 |
1,9E-05 |
-3,7E-05 |
6,6E-05 |
-0,00011 |
|
0,1 |
1,060831 |
0,078006 |
-0,00353 |
-0,00018 |
1,6E-05 |
1,1E-05 |
-1,8E-05 |
2,9E-05 |
-4,4E-05 |
|
|
0,2 |
1,138837 |
0,074475 |
-0,00371 |
-0,00016 |
2,7E-05 |
-7E-06 |
1,1E-05 |
-1,5E-05 |
|
|
|
0,3 |
1,213312 |
0,070764 |
-0,00388 |
-0,00014 |
2E-05 |
4E-06 |
-4E-06 |
|
|
|
|
0,4 |
1,284076 |
0,066889 |
-0,00401 |
-0,00012 |
2,4E-05 |
4,44E-16 |
|
|
|
|
|
0,5 |
1,350965 |
0,062877 |
-0,00413 |
-9,3E-05 |
2,4E-05 |
|
|
|
|
|
|
0,6 |
1,413842 |
0,058748 |
-0,00422 |
-6,9E-05 |
|
|
|
|
|
|
|
0,7 |
1,47259 |
0,054526 |
-0,00429 |
|
|
|
|
|
|
|
|
0,8 |
1,527116 |
0,050235 |
|
|
|
|
|
|
|
|
|
0,9 |
1,577351 |
|
|
|
|
|
|
|
|
|
Т.к. точность вычисления по условию нам не задана, интерполируем функцию полиномами 1 (линейная интерполяция), 5 и 8 порядков.
Ближайший к точке x=a=0.104 узел слева – x=0.1, поэтому полагаем x0=0.1. Точка a находится вначале таблицы, поэтому применяем «интерполяцию вперед», т.е. используем формулу (8), а для оценки погрешности формулу (12).
q= =0.04
P1(0.104)=y0+q ?y0
R1(0.104)
P5(0.104)=y0+q ?y0+ ?2y0+ ?3y0+ ?4y0+ ?5y0
R5(0.104)
P8(0.104)=y0+q ?y0+ ?2y0+ ?3y0+ ?4y0+ ?5y0+ + ?6y0+ ?7y0+
+ ?8y0
R8(0.104) .
Результаты интерполяции функции полиномом Ньютона методом «интерполяция вперед» представлены в таблице 3. Общий график функции и график функции в окрестности искомой точки представлены на рисунке 5 и рисунке 6, соответственно.
Как видно из графика (рис. 5) функция близка к линейной и линейная (n=2) интерполяция полиномом Ньютона даёт результат с достаточно высокой точностью. Что подтверждается графиком функции в окрестности точки х=0.104 (рис.6), все три графика практически совпадают.
Таблица 3
|
Х |
n |
Pn(x) |
E |
|
0.104 |
2 |
1.06395 |
6.7*10—5 |
|
0.104 |
5 |
1.06401 |
1.2*10—7 |
|
0.104 |
8 |
1.06402 |
1.9*10—7 |
Алгоритм программы интерполирования функции полиномом Ньютона представлен на рис. 7 и рис. 8. Программа включает в себя функцию «newton» принимающую массив табличных значений аргументов и функции, количество узлов интерполяции, значение точки «х», требуемую точность интерполирования результат вычисления полинома записывается в нулевой элемент двумерного массива.
|
Рис. 5. График заданной функции. |
|
Х |
|
Y |
|
Рис. 6. График функции, интерполированной полиномом Ньютона различного порядка, в окрестности точки х=0.104 |
|
Х |
|
Y |
|
Рис.7. Блок-схема алгоритма интерполирования «вперед» функции полиномом Ньютона |
|
newton |
|
i=0,N |
|
i=1,n |
|
Вывод arg(i,1), arg(i,2) |
|
Конец |
|
Начало |
|
Ввод N, X, E, h |
|
Ввод arg(i,2) |
|
arg(i,1)=i*h |
|
Вывод x, arg(0,0), arg(1,0) |
|
Табличные значения аргументов функции (х(i)) сохраняются в первом столбце (arg(i,1)) двумерного массива. Значения функции (y(i)) записываются во второй столбец (arg(i,2)) массива. Результат интерполирования и результат оценки погрешности сохраняются в нулевом (arg(0,0)) элементе и первой строке нулевого столбца (arg(1,0)) массива соответственно. |
|
Рис.8 Блок-схема алгоритма функции newton |
|
x=arg(i,1) |
|
newton(arg,x,n,e) |
|
i=0,n |
|
dX=x-arg(i,1) |
|
h=arg(2,1)-arg(1,1) |
|
(dX<h)&(dX>0) |
|
k=flag=0 |
|
flag=i |
|
k=i |
|
j=3,n+2 |
|
s=arg(flag,2) |
|
конец |
|
flag=0 |
|
j=k,n-j+2 |
|
arg(i,j)=arg(i+1,j-1)-arg(i,j-1) |
|
q=arg(k,1)/h |
|
m=1,n-k |
|
arg(0,0)=arg(k,2) |
|
i=1,j-2 |
|
s=s*((q-i+1)/i) |
|
arg(0,0)=arg(0,0)+s |
|
arg(max,max)=m; arg(1,0)=s |
|
| s |<e |
Текст программы вычисления многочлена Ньютона
// method_N.cpp : main project file.
#include «stdafx.h»
#include<iostream>
#include<conio.h>
#include<cmath>
#include<iomanip>
#define max 20
using namespace System;
using namespace std;
void newton (double arg[][max],double x,int n, double e) /*объявление функции «ньютон», входные данные:двумерный массив значений «Х» и «Y», значение «Х», введенное с клавиатуры, количество узлов интерполяции, заданная погрешность*/
double s,dX,q,h; //локальные переменные
int k, flag;
k=0;
flag=0;
h=arg[2][1]-arg[1][1];
for(int i=0;i<n;i++)
if(x==arg[i][1]) /*проверка равенства введенного значения аргумента «х» со значением из массива табличных аргументов «х[]» если равенство выполняется, в переменную flag записывается номер элемента массива табличных значений, равного введенному с клавиатуры «х»*/
flag=i;
dX=x-arg[i][1];
if((dX<h)&&(dX>0)) //нахождение ближайшего узла интерполяции
k=i;
if (flag==0) //если совпадений, введенного с клавиатуры значения, «х» с массивом табличных значений нет, вычисляются конечные разности
for(int j=3;j<=n+2;j++)
for(int i=k;i<=n-j+2;i++)
arg[i][j]=arg[i+1][j-1]-arg[i][j-1]; /* вычисление конечных разностей. Конечные разности записываются в двумерный массив
начиная с третьего столбца*/
q=(x-arg[k][1])/h;
arg[0][0]=arg[1][0]=arg[max][max]=0;
for(int m=1;m<=n-k;m++)
arg[0][0]=arg[k][2]; //инициализация первого слагаемого (y0) полинома
for(int j=3;j<=m+2;j++) //вычисление слагаемых полинома
{
s=arg[k][j]; //переменной s присваивается значение dY порядка j
for(int i=1;i<=j-2;i++)
s=s*((q-i+1)/i);
arg[0][0]=arg[0][0]+s; //вычисление значения функции в заданной точке
arg[1][0]=s; //в элементе первой строки нулевого столбца записывается приблизительная погрешность интерполирования
arg[max][max]=m; //в элемент последнего столбца нижней строки сохраняется порядок полинома
if (fabs(s)<e) //если заданная точность интерполирования достигнута процесс вычисления завершается
break;
else
s=arg[flag][2];
} //конец функции «ньютон»
int main() //начало основной программы
double arg[max][max],x,e;
float h;
int n;
cout<<«Enter quantity of figures N=»;
cin>>n; //ввод количества узлов интерполяции
cout<<endl;
cout<<«Enter value of argument X=»;
cin>>x; //ввод значения «х»
cout<<endl;
cout<<«Enter an interpolation error E=»;
cin>>e; //ввод заданной погрешности интерполяции
cout<<endl;
cout<<«Enter distance between knots h=»;
cin>>h; //ввод расстояния между узлами
cout<<endl;
for (int i=0;i<=n;i++)
{
arg[i][1]=i*h; //заполнение массива табличных значений аргумента функции «x», аргументы функции заполняют первый столбец двумерного массива
cout<<«Enter values of function Y»<<i<<«= «;
cin>>arg[i][2]; //заполнение массива табличных значений функции «y», значения функции заполняют второй столбец двумерного массива
cout<<«\n»;
newton(arg,x,n,e);
cout<<setw(5)<<«X(i)»<<setw(15)<<«Y(i)»<<endl;
for(int i=0;i<=n;i++)
cout<<setw(5)<<arg[i][1]<<setw(16)<<arg[i][2]<<«\n»; //вывод на экран табличных значений аргументов и функции
cout<<endl;
cout<<«X=»<<x<<» P(«<<arg[max][max]<<«)=»<<arg[0][0]<<» E=»<<arg[1][0]<<endl; // вывод результатов вычислений на экран
getch();
return 0;
2.2.2 Интерполяция «назад».
Если искомая точка находится в конце таблицы заданных значений, к примеру, a=x=0.865, то для того, чтобы использовать большее число узлов применяется интерполяция «назад», в которой используется вторая интерполяционная формула Ньютона (11).
Ближайший справа к точке a=x=0.865 узел справа х=0.9, поэтому полагаем хn=0.9. Используя интерполяционную формулу (11) и таблицу 2 составим полином Ньютона взяв два узла интерполяции. Для оценки погрешности используем формулу (13).
Для вычисления y(0.865) примем хn=0.9, уn=1.577351, тогда
q= = -0.35
P(0.865)=yn+q ?yn-1+ ?2yn-2 =
= 1.577351+(-0.35)* 0,050235+ -0,00429)= 1.5605739
R(0.865) 6.9*10-5=4.3*10-6
Выводы:
1. В формулах Ньютона в случае добавления узла все найденные члены сохраняются и появляется новое слагаемое, представляющее собой ни что иное, как поправку к уже вычисленному значению.
2. При интерполяции на малых участках слагаемые в формулах (10) и (11) будут расположены в порядке их малости, что облегчает использование формул Ньютона в вычислениях и позволяет судить о точности интерполяции.
3. Степень интерполирующего полинома существенно зависит от шага таблицы (чем меньше шаг, тем график функции более приближен к линейному, что позволяет использовать линейную интерполяцию)
4. Ограниченность применения формул Ньютона связанна с их пригодностью лишь для равноотстоящих узлов.
3. Интерполяция функций методом наименьших квадратов.
3.1. Теоретические основы метода
Областью применения является интерполяция экспериментальных данных, полученных со значительными погрешностями.
Наиболее распространен способ выбора функции ?(x)
?(x)=C0?0(x)+ C1?1(x)+…+ Cm?m(x) (14)
?0(x), ?1(x),… ?m(x) – базисные функции m?n,
=(C0,C1,…Cm) – вектор коэффициентов, определяемых из условия минимизации
где (n+1) — количество узловых точек.
Если погрешность исходных данных Е, то количество базисных функций выбирается так, чтобы ?Е.
Условие минимума функции S приводит к системе линейных уравнений относительно параметров C0,C1,…Cm. Система называется системой нормальных уравнений, ее матрица – матрицей Грама. Элементы матрицы Грама являются скалярными произведениями базисных функций
Для степенного базиса имеем систему уравнений:
столбец правых частей системы нормальных уравнений:
Возникающая при интерполяции ошибка (невязка) является наименьшим квадратичным отклонением:
где (n+1) – количество узлов, m – степень аппроксимирующего полинома, n+1?m.
Система линейных уравнений для нахождения коэффициентов многочлена ?0(x)=a0+a1x (линейная аппроксимация):
Введем обозначения – средние значения исходных данных.
Решением системы уравнений являются:
3.2. Составление аппроксимирующего полинома для выбранной функции и оценка точности аппроксимации.
Составим два варианта аппроксимации следующей функции для пяти узлов (n=5), используя метод наименьших квадратов:
|
Х |
Y |
|
0.0 |
0.979498 |
|
0.1 |
1.060831 |
|
0.2 |
1.138837 |
|
0.3 |
1.213312 |
|
0.4 |
1.284076 |
|
0.5 |
1.350965 |
|
0.6 |
1.413842 |
|
0.7 |
1.472590 |
|
0.8 |
1.527116 |
|
0.9 |
1.577351 |
3.2.1. Линейная аппроксимация.
Найдем средние значения исходных данных:
= 0.2;
= 1.1353108
Далее найдем коэффициенты многочлена:
= 0.076367/0.1 = 0.76367
= 1.1353108-0.76367*0.2= 0.9825768.
Итак, полином для линейной аппроксимации равен:
?(x)= 0.76367+0.9825768*x
|
x |
f(x) |
?(x) |
|
0.0 |
0.979498 |
0,76367 |
|
0.1 |
1.060831 |
0,861928 |
|
0.2 |
1.138837 |
0,960185 |
|
0.3 |
1.213312 |
1,058443 |
|
0.4 |
1.284076 |
1,156701 |
Используя формулу(18), вычислим невязку:
P= = 0.1779157
На рис. 9 график исходной функций и функции, найденной линейной аппроксимацией.
Как видно из результатов вычисления невязки и графика функций (рис.9) видно, что найденная функция не соответствует исходной. Качество аппроксимации неудовлетворительное, следовательно, для заданной функции линейная аппроксимация методом наименьших квадратов неприменима.
|
Рис. 9. график исходной функций и функции, найденной линейной аппроксимацией.
|
|
Х |
|
Y |
3.2.2. Квадратичная аппроксимация.
Запишем в таблицу 4 элементы матрицы Грама и столбец свободных членов.
Таблица 4.
|
|
x0 |
x1 |
x2 |
x3 |
x4 |
y |
xy |
x2y |
|
1 |
0 |
0 |
0 |
0 |
0,979498 |
0 |
0 |
|
|
1 |
0,1 |
0,01 |
0,001 |
0,0001 |
1,060831 |
0,106083 |
0,010608 |
|
|
1 |
0,2 |
0,04 |
0,008 |
0,0016 |
1,138837 |
0,227767 |
0,045553 |
|
|
1 |
0,3 |
0,09 |
0,027 |
0,0081 |
1,213312 |
0,363994 |
0,109198 |
|
|
1 |
0,4 |
0,16 |
0,064 |
0,0256 |
1,284076 |
0,51363 |
0,205452 |
|
|
? |
5 |
1 |
0,3 |
0,1 |
0,0354 |
5,676554 |
1,211475 |
0,370812 |
Построим матрицу Грамма:
Используя метод Крамера решим систему линейных уравнений:
= 0,0007
= 0,000685621
= 0,000582484
= -0,000123345
С0 = 0,9794592
С1 = 0,8321198
С2 = -0,1762071
Окончательно искомая аппроксимирующая функция примет вид:
?(х)=0,9794592+0,8321198*х-0,1762071*х2.
Результаты аппроксимации функции сведены в таблицу 5.
Таблица 5.
|
x |
y |
= ?(х) |
|
0 |
0,979498 |
0,979459 |
|
0,1 |
1,060831 |
1,060909 |
|
0,2 |
1,138837 |
1,138835 |
|
0,3 |
1,213312 |
1,213237 |
|
0,4 |
1,284076 |
1,284114 |
Вычислим невязку:
P= = 5.432*10-5
На рисунке 10 представлены графики исходной и аппроксимирующей функций. Блок-схема алгоритма реализации аппроксимации функции методом наименьших квадратов представлена на рис.10. Основная программа включает в себя функции «gram» (рис.11) — строит матрицу Грамма и «matrix» (рис.12) — реализующую вычисление систем линейных уравнений методом Холесского. Функция «matrix» в свою очередь, включает в себя функцию «holess» (рис.14) — преобразующую массив исходных значений (исходную матрицу) в массив новых коэффициентов. И функцию «res» (рис.13) на основании массива, преобразованного функцией «holes», вычисляющую и заполняющую массив готовых решений. Все функции имеют не возвращаемый тип void, т.к. работают с исходным массивом значений.
|
Рис. 10. график исходной функций и функции, найденной квадратичной аппроксимацией.
|
|
Х |
|
Y |
|
Конец |
|
s=s+(y1[i]-y[i][0])2 |
|
вывод x(i,1), y(i,0),y1(i) |
|
j= |
|
Начало |
|
ввод n |
|
i= |
|
ввод x(i,1) |
|
ввод y(i,0) |
|
gram |
|
matrix |
|
i= |
|
s=0 |
|
x(i,j)=x(i,1)j |
|
j= |
|
y(i,j)=y(i,0)*x(i,1) |
|
j= |
|
s=s+c[j]*x[i][j-1] |
|
j= |
|
y1(i)=s |
|
s=0 |
|
p= |
|
вывод p |
|
Рис.10. Блок-схема алгоритма метода наименьших квадратов |
|
j= |
|
x(n+1,j)=0; y(n+1,j)=0 |
|
gram (x(i,j),y(i,j),n,a(i,j)) |
|
i= |
|
x(n+1,j)=x(n+1,j)+x(i,j) |
|
j= |
|
y(n+1,j)=y(n+1,j)+y(i,j) |
|
i= |
|
k=0 |
|
i= |
|
a(i,j)=x(n+1,j-1+k) |
|
j= |
|
a(i,n-1)=y(n+1,i-1) |
|
k=k+1 |
|
Конец |
|
Рис.11. Блок-схема алгоритма функции «gram» |
|
matrix (a(i,j),n,c(i)) |
|
i= |
|
c(i)=0 |
|
holess |
|
res |
|
Конец |
|
Рис.12. Блок-схема алгоритма функции «matrix» |
|
res ((arg(i,j),n,x(i) |
|
i=
|
|
n=1 |
|
s=0 |
|
k=
|
|
s=s+arg(i,k)*x(k) |
|
s=s+arg(i,k)*x(k) |
|
Конец |
|
Рис. 13. Блок-схема алгоритма функции «res» |
|
Конец |
|
||arg(t,m)|?max |
|
t= |
|
max=|arg(m,m)| |
|
holes(arg(i,j),n) |
|
m= |
|
m=1 |
|
j= |
|
arg(1,j)=arg(1,j)/arg(1,1) |
|
m ? 2 |
|
max=|arg(t,m)|; i1=t |
|
i= |
|
a1(i)= arg(m,i) |
|
i= |
|
arg (m,i)= arg(i1,i) |
|
arg (i1,i)= a1(i) |
|
i= |
|
p=0 |
|
k= |
|
p=p+arg(i,k)*arg(k,m) |
|
arg(i,m)=arg(i,m)-p |
|
j= |
|
s=0 |
|
k= |
|
s=s+arg(m,k)*arg(k,j) |
|
arg(m,j)=(arg(m,j)-s)/arg(m,m) |
|
Рис.14. Блок-схема алгоритма функции «holess» |
Текст программы реализующей аппроксимацию методом наименьших квадратов
// min_square.cpp : main project file.
#include «stdafx.h»
#include<iostream>
#include<conio.h>
#include<cmath>
#include<iomanip>
#define size 20
using namespace System;
using namespace std;
//объявление прототипов функций
void matrix (double[][size],int,double[]); /*вычисляет значение матрицы методом Холесского, включает в себя две функции: holess, и res*/
void holess(double[][size],int);
void res(double[][size], int,double[]);
//функция построения матрицы Грама
void gram(double[][size],double[][size],int, double[][size]);
//начало основной программы
int main()
{
/*объявленные переменные: в массивы «х» и «у» накапливаются табличные значения аргументов и функции. В массив «с» сохраняются значения коэффициентов (результат решения матрицы Грама). Массив «а» хранит матрицу Грама, «у1»-значения функции, полученные путем вычисления аппроксимирующего многочлена. s-хранит промежуточные результаты вычисления функции апроксимирующим многочленом. р-невязка.n-количество табличных значений аргументов*/
double x[size][size],y[size][size],c[size],a[size][size],y1[size],s,p;
int n;
cout<<«Enter quantity of knots of interpolation N=»; //ввод n
cin>>n;
cout<<endl;
for (int i=1;i<=n;i++) //начало цикла ввода табличных значений
cout<<«Inter x(«<<i<<«)=»;
cin>>x[i][1]; /*табличные значения аргумента сохраняются в первом
столбце многомерного массива «х»*/
for(int j=0;j<=(n-1);j++)
x[i][j]=pow(x[i][1],j); /*табличное значение «х» первого столбца i-й
строки возводится в j-ю степень и сохраняется в j-м столбце*/
cout<<«Inter y(«<<i<<«)=»;
cin>>y[i][0]; /*сохранение табличных значений функции в нулевом
столбце многомерного массива «у» */
for(int j=0;j<=(n-3);j++)
y[i][j]=y[i][0]*x[i][j]; /*в j-й столбец i-й строки массива «у» записывается произведение значения 0-го столбца массива «у» (табличное значение у) и х^j */
cout<<endl;
gram(x,y,n,a); //заполнение матрицы Грама
matrix(a,n-2,c); //вычисление матрицы методом Холесского
for(int i=1;i<=n;i++)
s=0;
for(int j=1;j<=n-2;j++)
s=s+c[j]*x[i][j-1]; //вычисление функции аппроксимирующим многочленом
y1[i]=s; //запись вычисленного значения функции в массив выходных данных
s=0;
cout<<setw(2)<<«X»<<setw(12)<<«Y»<<setw(19)<<«F»<<endl;
for(int i=1;i<=n;i++)
s=s+pow((y1[i]-y[i][0]),2); //нахождение разницы между табличным и вычисленным значением функции
cout<<setw(4)<<x[i][1]<<setw(14)<<y[i][0]<<setw(19)<<y1[i]<<endl; //вывод табличных и вычисленных значений
p=sqrt(s/(n+1)); //вычисление невязки
cout<<endl;
cout<<setw(2)<<«P=»<<p<<endl; //вывод невязки на экран
getch();
return 0; //конец функции main
/*начало функции, строящей матрицу Грама. Входные данные: двумерный массив табличных, а также возведенных в степень n-1 значений аргументов «х»; двумерный массив табличных значений, а так же результатов произведений аргумента и функции «у»; количество узловых точек. Выходные данные двумерный массив «а» — матрица Грама*/
void gram(double x [][size],double y[][size],int n, double a[][size])
int k;
for(int j=0;j<=(n-1);j++)
x[n+1][j]=y[n+1][j]=0;
for(int i=1;i<=n;i++)
x[n+1][j]=x[n+1][j]+x[i][j]; //вычисление значений членов левой части
матрицы Грама
cout<<endl;
for(int j=0;j<=n-3;j++)
for(int i=1;i<=n;i++)
y[n+1][j]=y[n+1][j]+y[i][j]; //вычисление значений членов правой части матрицы Грама
k=0;
for(int i=1;i<=n-2;i++)
for(int j=1;j<=(n-2);j++)
a[i][j]=x[n+1][j-1+k]; //заполнение массива «а» значениями аргументов
a[i][n-1]=y[n+1][i-1]; //заполнение столбца правых частей системы
нормальных уравнений
k++;
void matrix (double a[][size],int n,double c[])
for(int i=1;i<=n;i++)
c[i]=0;
holess(a,n);
res(a,n,c);
//функция реализует метод Холесского с частичной перестановкой для выбора главного элемента
//Входные данные массив коэффициентов «аrg», количество уравнений системы n
//Выходные данные массив новых коэффициентов «аrg».
void holess(double arg[][size],int n)
double max,p,s,a1[size];
int i1;
for(int m=1;m<=n;m++)
{
if(m==1)
for(int j=2;j<=(n+1);j++)
arg[1][j]=arg[1][j]/arg[1][1];
if(m>=2)
max=fabs(arg[m][m]);
for(int t=m;t<=n;t++)
if(fabs(arg[t][m])>=max) //выбор главного элемента
{
max=fabs(arg[t][m]);
i1=t;
}
for(int i=1;i<=(n+1);i++)
a1[i]=arg[m][i]; //a1(n+1) — массив коэффициентов
for(int i=1;i<=(n+1);i++) //замена ведущей строки на строку с выбранным главным элементом
arg[m][i]=arg[i1][i];
arg[i1][i]=a1[i];
for(int i=m;i<=n;i++)
p=0;
for(int k=1;k<=(m-1);k++)
p=p+arg[i][k]*arg[k][m];
arg[i][m]=arg[i][m]-p; //формирование очередного a(i,m) столбца
for(int j=(m+1);j<=(n+1);j++)
s=0;
for(int k=1;k<=(m-1);k++)
s=s+arg[m][k]*arg[k][j];
arg[m][j]=(arg[m][j]-s)/arg[m][m]; //формирование очередной a(m,j)строки
}
//функция формирования решения системы уравнений
//входные данные: двумерный массив коэффициентов «arg», n-число уравнений
//выходные данные:массив решений системы х(i)
void res(double arg[][size], int n,double x[])
double s;
for(int i=n;i>=1;i—)
if(n==1)
x[i]=arg[n][n+1];
s=0;
for(int k=i+1;k<=n;k++)
s=s+arg[i][k]*x[k];
x[i]=arg[i][n+1]-s;
} //конец цикла
} //конец функции res
Выводы:
Результаты расчетов, а также построенные графики позволяют сделать вывод, что квадратичная аппроксимация в данном случае предпочтительнее, т.к. при квадратичной аппроксимации график аппроксимирующей функции на рассчитанном участке с достаточно высокой точностью повторяет график исходной функции. Невязка при квадратичной аппроксимации значительно меньше, чем при линейной.
Заключение:
В настоящей работе были рассмотрены три способа аппроксимации одной и той же функции. У всех трех способов есть своя область применения, достоинства и недостатки. Применимо к исследуемой функции, учитывая равный интервал между значениями аргументов, малый шаг, а так же то, что функция близка к линейной, наиболее предпочтительным для ее аппроксимации является метод Ньютона. Т.к. уже линейная интерполяция дала результат с высокой точностью.
Список литературы:
- Гловацкая А.П. Методы и алгоритмы вычислительной математики. Уч.пос. для ВУЗОВ.-М. Радио и связь, 1999.-408с.
- Гловацкая А.П. Сборник задач для курсовой работы по курсу Информатика. МТУСИ. 2006, 32с.
- Дейтел Х.М. Как программировать на С++: Пятое издание. — С.Пб. ООО «Бином-Пресс», 2011.-1456с.
- Пахомов Б.И. С/С++ и MS Visual C++ 2010 для начинающих – СПб.: БХВ — Петербург, 2011. – 736 с.
