Ращиков В.И

Задание № 3

ИНТЕРПОЛИРОВАНИЕ

Цель работы — изучение методов интерполирования, построение интерполяционного многочлена Ньютона.

ПОСТАНОВКА ЗАДАЧИ

Численное моделирование большинства физических задач сопряжено, как правило, с необходимостью учета факторов, которые не могут быть описаны аналитически. Имеется лишь ряд экспериментальных зависимостей, полученных в фиксированном числе точек интересующего нас диапазона переменных. Так, при решении широко распространенной задачи о динамике макро- и микрообъектов во внешних гравитационных или электромагнитных полях информацию о поле часто бывает невозможно получить в виде аналитических функций без ввода дополнительных упрощающих предположений, которые могут существенно повлиять на результат. В этом случае необходимо прибегнуть к экспериментальным характеристикам, причем эксперимент может быть проведен лишь конечное число раз. Таким образом, мы приходим к физической задаче, в которой ряд функций задан на конечном числе точек xi фиксированной области изменения аргумента xi [a, b]. Численный метод, однако, может потребовать знания этих функций для всех значений аргумента из этой области. В этом случае возникает задача восстановления функции y(x) для всех значений x [a, b], если известны ее значения в некотором фиксированном числе точек xi этого отрезка.

Наиболее простым и распространенным способом решения этой задачи является интерполяция между соседними значениями, которая сводится к построению функции ϕ(х), совпадающей с функцией y(x) в точках xi, т. е.

ϕ(xi) = y(xi) = yi, i = 0, 1, 2, …, n,

где n + 1 — число заданных на отрезке [a, b] точек, а xi — узлы интерполяции.

21

При выборе интерполирующей функции ϕ(x) необходимо ограничить поиск функциями, легко и быстро вычисляющимися на компьютере, так как их, как правило, приходится вычислять многократно.

Существует много интерполяционных многочленов и способов их построения, пригодных для различного расположения узлов. При построении интерполяционных многочленов обычно подразумевается, что множество используемых узлов известно. Однако часто бывает известна лишь требуемая точность, а число узлов не фиксировано.

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

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

Задание включает следующие этапы.

1.

Вычислить таблицу значений y j

= y(xj ) заданной функции

у(х) в равноотстоящих узлах

xj = jh

( j = 0,1,...,n),

h = 1/ n,

отрезка [0,1].

2.

Составить таблицу первых разностей функции

y(x j+1, x j ) =

y j+1 − y j

=

y j+1 − y j

( j = 0,1,..., n −1).

h

3.

x j+1 − x j

Составить таблицу вторых разделённых разностей

y(x

j+2

, x

j+1

, x

j

) = y(x j+2 , x j+1 ) − y(x j+1, x j ) = y j+2 −2 y j+1 + y j

x j+2 − x j

2h2

( j = 0,1,...,n −2).

4.

По этим таблицам, используя интерполяционный многочлен

Ньютона второго порядка

P(x) = yj +(x − xj ) y(xj+1, xj ) +(x − xj )(x − xj +1 ) y(xj+2 , xj+1, xj ) =

= y

j

+(x − x

j

)

y j+1 − y j

+(x − x

j

)(x − x

j+1

)

y j+2 −2yj +1 + y j

,

h

2h2

22

вычислить значения Р(х) в точках (узлах) с полуцелыми индексами xj+1/2 = ( j +1/ 2)h ( j = 0,1,...,n −2).

5. Найти погрешность интерполирования в этих узлах

εj +1/2 = y(xj+1/ 2 ) − P(xj+1/ 2 ) ( j = 0,1,..., n −2),

максимальную погрешность εmах, средний квадрат погрешности и среднеквадратичную погрешность εm:

							1	n−2
	εmax = max	εj+1/2	,	ε2 =				∑ε2j+1/ 2 , εm = ε2 .
	j						n −1 j=0

6. Исследовать, как изменяются погрешности εmах и εm с изменением n.

ВАРИАНТЫ ФУНКЦИЙ у(х)

	(a,b =1÷5; k, m =1÷4;		n = 20 ÷100)
1) y(x) =sink (πxm );		2) y(x) = cosk (πxm );
3) y(x) = sink (πx1/ m );		4) y(x) = cosk (πx1/ m );
5) y(x) = tgk (πxm / 4);		6) y(x) = tgk (πx1/ m / 4);
7) y(x) = ±ax4 ±bx3 ;		8) y(x) = (a + bxm )−k ;
9)	y(x) = (a +bxm )−1/ k ;	10)	y(x) = (a +bx1/ m )−k ;
11)	y(x) = (a +bx1/ m )−1/ k ;	12)	y(x) = xk /(a +bxm );
13) y(x) = xk /(a +bx)m ;		14)	y(x) = x1/ k /(a +bx)m ;
15) y(x) = xk /(a +bx)1/ m ;		16)	y(x) = x1/ k /(a +bx)1/ m ;
17) y(x) = (a + x2 )k / (b + x4 )m ;		18) y(x) = (a + x2 )1/k / (b + x4 )m ;
19) y(x) = (a + x2 )k / (b + x4 )1/ m ;		20) y(x) = (a + x2 )1/ k / (b + x4 )1/m ;
21) y(x) = xk / (a +bxm );		22) y(x) = x1/k / (a +bxm );
23) y(x) = xk / (a +bx1/ m );		24) y(x) = x1/ k / (a +bx1/m );
25) y(x) = lnk (1 + xm );		26) y(x) = ln1/k (1 + xm );

23

27)y(x) = lnk (1+ x1/m );

29)y(x) = xk e−2 x2 ;

31)y(x) = xk e−2 x1/m ;

33)y(x) = e−8( x−0.5)2 ;

35)y(x) = x1/k e−8( x−0.5)2 ;

37) y(x) = (a + bx1/ m )1/k ; 39) y(x) = (a + bxm )−1/k ;

41) y(x) = arcsin[cosk (0.5πxm )]; 43) y(x) = arcsin[cosk (0.5πx1/m )]; 45) y(x) = arctg(a ± bxm );

47) y(x) = sh(a ± bx)−m ; 47) y(x) = ch(a ± bx)−m ;

28)y(x) = ln1/k (1+ x1/m );

30)y(x) = x1/k e−2 x2 ;

32)y(x) = x1/k e−2 x1/m ;

34)y(x) = xk e−8( x−0.5)2 ;

36)y(x) = (a + bxm )1/k ;

38)y(x) = (a + bxm )−k ;

40)y(x) = (a + bx1/m )−1/k ;

42)y(x) = arccos[sink (0.5πxm )];

44)y(x) = arccos[sink (0.5πxm )];

46)y(x) = arctg(a ± bx1/m );

48)y(x) = sh(a ± bx)−1/m ;

50)y(x) = ch(a ± bx)−1/m .

ПРОГРАММИРОВАНИЕ

Блок-схема программы представлена на рис. 3.1. Основу программы составляют три последовательных цикла — блоки 3-4- 5, 6- 7-8, 9-10-11. Для хранения вычисляемых в этих циклах значений функции, первых и вторых разностей, а также погрешностей следует описать соответствующие массивы.

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

y(x) = x2 ,

для которой εmax и εm должны обращаться в ноль.

24

1

2

3

4

5

6

7

8

Начало

Начальные

данные

Цикл по j =1,n

yj

Цикл по j

Цикл по j =1,n-1

y(xj+1, xj)

Цикл по j

		Цикл
9
		по j =1,n-2

y( x j+ 2 , x j+1 , x j )

10P( x j +1/ 2 ),ε j +1/ 2 ,

εmax , ∑ε2j +1/ 2

11Цикл по j

12ε2 ,εm

Графики 13 Результаты

14 Конец

Рис.3.1. Блок-схема программы интерполяции

25

СОДЕРЖАНИЕ ОТЧЕТА

Отчет должен содержать:

•формулу и график функции у(х)для конкретного варианта;

•текст программы;

• таблицу погрешностейεj +1/ 2 ( j = 0,1,....., n − 2);

• значения εmax и εт.

КОНТРОЛЬНЫЕ ВОПРОСЫ

1.Как ставится задача интерполирования?

2.Какие интерполяционные многочлены вы знаете?

3.Как определяются разделенные разности различных поряд-

ков?

4. Как строился интерполяционный многочлен Ньютона?

5.Какова погрешность (остаточный член) интерполяционного многочлена?

6.Как можно практически оценить погрешность интерполирования?

26

Задание № 4

АППРОКСИМАЦИЯ

Цель работы: изучение аппроксимации функций на примере метода наименьших квадратов.

ПОСТАНОВКА ЗАДАЧИ

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

Рис. 4.1. Иллюстрация аппроксимации и интерполяции

Будем аппроксимировать функции многочленом степени m:

ϕ(x) = c0 + c1x + c2x2 + … + cmxm,

коэффициенты которого ci подберем так, чтобы минимизировать отклонение многочлена от данной функции.

27

Воспользуемся среднеквадратичным приближением функции y(x) многочленом ϕ(x) на множестве (xi, yi), (i = 0, 1,2, …, n), при котором мерой отклонения является величина S, равная сумме квадратов разностей между значениями многочлена и функции в данных точках:

n	n
S = ∑εi2	= ∑[ϕ(xi , c0 , c1, ..., cm ) − yi ]2 .
i=0	i=0

Для построения аппроксимирующего многочлена нужно подобрать коэффициенты c0, c1, …, cm так, чтобы величина S была наименьшей. В этом и состоит метод наименьших квадратов. Если отклонение εi подчиняется нормальному закону распределения, то полученные таким образом значения параметров наиболее вероятны. Как уже упоминалось выше, среднеквадратичное приближение сглаживает неточности функции, давая правильное представление о ней.

Поскольку ci выступают в роли независимых переменных функции S, то минимум найдем, приравнивая нулю частные производные по этим переменным:

	∂S	= 0,	∂S	= 0,	∂S	= 0, ...,	∂S	= 0,
	∂c		∂c		∂c		∂c
	0		1		2		n

т. е. приходим к системе уравнений для определения сi. Если в качестве аппроксимирующей функции взять многочлен, то выражение для квадратов отклонений примет вид:

n

S = ∑(c0 +c1xi +c2 xi2 +... +cm xim − yi )2 .

i=0

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

	∂S		n
			= 2∑(c0 + c1xi	+... + cm xim − yi ) = 0;
	∂c0
			i=0
	∂S		n
			= 2∑(c0 + c1xi	+... +cm xim − yi )xi = 0;
	∂c1
			i=0
	.......
	∂S		n
			= 2∑(c0 + c1xi	+... +cm xim − yi )xim = 0.
	∂cm
			i=0

28

Собирая коэффициенты при неизвестных c0, c1, …, cm, получаем систему уравнений:

	n	n	n	n
(n + 1)c0 + c1 ∑xi + c2 ∑xi2 + ... + cm ∑xim =				∑ yi ;
	i=0	i=0	i=0	i=0
n	n	n	n	n
c0 ∑xi	+ c1 ∑xi2 + c2	∑xi3 + ... + cm ∑xim+1 =		∑xi yi ;
i=0	i=0	i=0	i=0	i=0
.........
n	n	n	n	n
c0 ∑xim + c1 ∑xim+1 + c2 ∑xim+2 + ... + cm ∑xi2m = ∑xim yi .
i=0	i=0	i=0	i=0	i=0

Решая систему, находим неизвестные параметры c0, более компактном виде можно записать:

b00c0 + b01c1 b10c0 + b11c1

bm0c0 + bm1c1

+… + b0mcm = a0;

+… + b1mcm = a1;

…

+… + bmmcm = am,

c1, …, cm. В

(4.1)

n

если ввести обозначения bkl = ∑xik +l

i=0

n

, ak = ∑xik yi ; k,l = 0,1, ..., m.

i=0

В данной работе, для того чтобы облегчить решение системы,

ограничимся значением		m =2.
Обозначим чертой усреднение по множеству узлов хi
					1		n
				=			∑ui
		u
и введем также обозначения:					n + 1 i=0
mk =			(k = 1, 2				...), Kii =		.
	xk							xi yi
Тогда систему (4.1) можно записать в виде:
c0 + m1c1 + m2c2 = K01;
m1c0 + m2c1 + m3c2 = K11;

m2c0 + m3c1 + m4c2 = K21.

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

29

		s11 =1, s12 = m1,s13 = m2 ;
		s		=	m	−m2		, s	= (m				−m m ) / s ;
		22			2	1		23				3		1	2		22
		s = m −(s2						+ s2		) = m −(m2						+ s2			);
		33			4	13		23					4		2		23
				z1 = K01, z2 = (K11 −m1K01 ) / s22 ;
				z3 =[K21 −(s13z1 + s23z2 )] / s33;
c	2	=	z3		,c =		z2 − z23c2			, c	0	= z		−(s c			+ s c				).
			s33		1		s22						1		12	1		13		2

(Для системы трех линейных уравнений вместо метода квадратных корней несложно воспользоваться известной схемой Крамера.)

Найдя коэффициенты c0, с1, с2, вычислим значения полинома ϕ(xj) (j = 0,1,..., n) и погрешность аппроксимации

εmax = max εj , εj = y(xj ) −φ(xj ), εm = ε2 ,

j

где средний квадрат погрешности ε2 = S / (n +1).

ВАРИАНТЫ ЗАДАНИЙ

Параметры: a = 0, b =1, n =10 ÷50.

Узлы: xj = jh( j = 0,1,.....,n), h =1 / n.

ВАРИАНТЫ ФУНКЦИЙ у(х)

1)y(x) = sin xq (q =1÷3);

2)y(x) = sin x1/q (q = 2,3);

3)y(x) = ex ;

4)y(x) = ln(1+ xq )(q =1÷3);

5)y(x) = cos x1/q (q = 2,3);

6)y(x) = cos xq (q =1÷3);

7)y(x) = ex1/ q (q = 2,3);

8)y(x) = ln(1+ x1/q )(q =1÷3);

30