Вычислит.матем_пособие

Pn ( x )

1 ha1

2

ha2 ( x

xn 1 )

... n

han ( x

xn 1)( x xn 2 )...( x x1 ).

Отсюда, полагая х=хn-1

получим a

yn 1

. Из выражения для второй

1

h

конечной

разности

найдем

а2:

a2

2 yn

2

.

Общая

формула для

2! h 2

коэффициента аi имеет вид ai

i yn

i

.

i! hi

Подставим

эти коэффициенты в

формулу многочлена и получим

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

P ( x )

y

yn 1

( x

x

) ...

n y0

( x

x

)...( x

x ).

n

n

n

n

h

n! hn

1

На практике используют формулу Ньютона в другом виде. Положим q=(x-xn)/h. Тогда

P ( x ) y

n

q y

n 1

q( q

1) 2

y

n 2

...

q( q

1)...( q

n 1)

n

y

0 .

n

2!

n!

5.3 Интерполирование сплайнами

Многочлен

Лагранжа

или

Ньютона

на

всем

отрезке

a,b с

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

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

Определение. Сплайн - функцией называют кусочно-полиномиальную функцию, определенную на отрезке a,b и имеющую на этом отрезке некоторое число непрерывных производных.

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

51

52

5.3.1 Построение кубического сплайна
Пусть на	отрезке		a,b	в	узлах		сетки		заданы значения		некоторой
функции f (x) , т.е. a x		0	x	x	2	... x	n	b ,	y	f (x ) (i= 0,1,…, n).
			1						i	i
Сплайном,	соответствующим этим узлам функции f (x)										называется

функция S(х), которая:

1) на каждом частичном отрезке является многочленом третьей степени;

2) функция S( x ) и ее первые две производные

S ( x ), S ( x )

непрерывны

на

a,b ;

3) S(xi )

f (xi ) .

На каждом

частичном

отрезке

xi 1 ,xi

будем

искать

сплайн

S(x)

Si (x) , где

Si (x)

многочлен третьей степени

Si (x)

ai

bi (x

xi )

ci

(x

xi

) 2

di

(x

xi ) 3 .

(5.8)

2

6

То есть для

x

xi 1 ,xi

нужно построить такую функцию Si ( x ) , где

ai ,bi ,ci ,di

подлежат определению.

Для всего отрезка интерполирования

a,b , таким образом, необходимо

определить

4 n

неизвестных

коэффициента.

S ( x ) b

c

( x x

)

di

( x x

)

2

,

i

i

i

i

2

S ( x ) ci

di ( x xi ),

Si ( x )

ai

yi .

Доопределим

a0

f (x0 ) y0 .

Требование непрерывности функции S(x)

приводит к условиям Si (xi )

Si

1 (xi ), (i=0, 1,…,n-1).

Отсюда из (5.8) получаем следующие уравнения:

ai

ai

bi 1( xi

xi 1 )

ci 1

( xi

xi

1 )

2

di 1

( xi

xi 1 )

3

(i= 1,2,…,n-1).

1

2

6

Введем шаг интерполирования hi

xi

xi

1 . Тогда последнее равенство

h2

h3

можно

переписать

в виде

hi

bi

i

ci

i

di

fi

fi

1 (i= 1,2,…,n). Из

2

6

h2

непрерывности первой

производной

следует

hi

ci

i

di

bi

bi

1 (i=

2

2,3,…,n), а из непрерывности второй производной

h d

i

c

c

1

(i=

i

i

i

2,3,…,n).

53

Объединив все три вида уравнений, получим систему из 3n-2 уравнений относительно 3n неизвестных bi ,ci ,di . Два недостающих уравнения получим, задав граничные условия для функции S(x). Для этого воспользуемся граничными условиями для сплайн-функции в виде

S ( a ) S ( b )

0 (концы гибкой линейки свободны).

Тогда получим систему уравнений

hi

di

ci

ci

1 ,c0

cn

0,( i

1,2,..., n )

h2

hi

ci

i

di

bi

bi 1 ,( i

2,3,..., n )

2

(5.9)

h2

h3

hi

bi

i

ci

i

di

fi

fi

1 ,( i

1,2,..., n ).

2

6

Решая систему методом подстановки (исключаем из (5.9) неизвестные bi,di), получим систему:

h c

i 1

2( h

h

) c

i

h c

6 (

yi 1

yi

yi yi

1

)

i

i

i 1

i 1 i 1

hi

hi

1

(5.10)

c0 cn

0

(i= 1,2,…,n-1).

Система (5.10) имеет трехдиагональную матрицу. Эта система может быть решена методом прогонки или Гаусса. После ее решения коэффициенты

сплайна di ,bi определим через коэффициенты сi с помощью явных формул

di	ci ci 1	,
	hi

	h		h2			y	i	y	i 1
b	i	c	i	d	i					(i= 1,2,…,n).
i	2	i	6					hi

5.3.2 Сходимость процесса интерполирования кубическими сплайнами

Доказывается, что при неограниченном увеличении числа узлов на

одном и том же отрезке		a,b	S(x)			f (x) . Оценка погрешности
интерполяции R( x ) f ( x ) S( x )			зависит			от	выбора сетки и степени
гладкости функции f(x).
При равномерной сетке xi		a	i h (i=0,1,…,n)
					M 4		h4
	f ( x )	Sh	( x )					,
						8

54

где M

4

max | f IV ( x ) |.

[ a ,b ]

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

1.

Если функция периодическая, то используется тригонометрическая

интерполяция

с периодом

l,

которая

строится

с

помощью

n

kx

kx

T ( x )

a

a

cos

b sin

тригонометрического многочлена

0

k

,

n

l

k

l

K 1

коэффициенты

которого находятся

из

системы

Tn ( xi )

f ( xi )

(i=

1,2,…,2n+1).

2. Выделяют приближение функций рациональными, дробно – рациональными и другими функциями. В данной книге эти вопросы не рассматриваются.

5.4 Аппроксимация функций методом наименьших квадратов

К такой задаче приходят при статистической обработке экспериментальных данных с помощью регрессионного анализа. Пусть в результате исследования некоторой величины x значениям x1 ,x2 ,x3 ,...,xn

поставлены в соответствие значения y1 , y2 , y3 ,..., yn некоторой величины у.

Требуется подобрать вид аппроксимирующей зависимости y=f(x), связывающей переменные х и у. Здесь могут иметь место следующие случаи. Во-первых: значения функции f(x) могут быть заданы в достаточно большом количестве узлов; во-вторых: значения таблично заданной функции отягощены погрешностями. Тогда проводить приближения функции с помощью многочлена нецелесообразно, т.к.

-это неудобно делать, поскольку число узлов велико и пришлось бы строить несколько интерполяционных многочленов;

-построив интерполяционные многочлены, мы повторили бы те же самые ошибки, которые присущи таблице.

Будем искать приближающую функцию из следующих соображений: 1) приближающая функция не проходит через узлы таблицы и не повторяет ошибки табличной функции; 2) чтобы сумма квадратов отклонений приближающей функции от таблично заданной была минимальной.

у	отклонения

55

х0

х1

хn-1 хn х

Рисунок 6 – Графическое изображение отклонений

Рассмотрим линейную задачу наименьших квадратов.

Пусть даны функции 0 ( x ), 1( x ),..., m ( x ), назовем их базисными

функциями. Будем искать приближающую (аппроксимирующую) функцию в виде линейной комбинации

y Фm( x ) c0 0( x ) c1 1( x ) ... cm m( x ) .

(5.11)

Такая аппроксимация называется линейной, а Фm(х) – обобщенный многочлен. Согласно критерию метода наименьших квадратов вычислим сумму квадратов отклонений таблично заданной функции от искомого многочлена в узлах:

	n		n
m	( yi	Фm ( xi ))2	( yi c0 0( xi ) ... cm m ( xi ))2 .	(5.12)
	i 0	i	0
Но нам неизвестна степень обобщенного многочлена. Подберем ее так,
чтобы	m было наименьшим и:

-аппроксимирующая кривая не проходила через узлы таблицы;

-получить приближение с заданной степенью точности.

Выражение m можно рассматривать как функцию от неизвестных

c0 ,..., cm . Нас интересует, при каких значениях c0 ,..., cm , значение m будет

минимально.

Для этого воспользуемся условием существования экстремума, а именно, найдем частные производные от m по всем переменным c0 ,..., cm и приравняем их к нулю. Получим систему вида:

56

			n
m	2		( yi	c0	0 ( xi )	...	cm	m ( xi ))	0 ( xi )	0
c0
		i	0
. . . . . . . . . . . . . . . . ..
			n							(5.13)
m	2		( yi	c0	0 ( xi )	...	cm	m ( xi ))	m ( xi )	0
cm
		i	0

Система (5.13) - система линейных уравнений относительно c0 ,..., cm .

Введем определение, чтобы лучше записать (5.13).

Определение. Скалярным произведением функции f на g на множестве

n

точек x0 ,..., xn

называется выражение ( f , g )

f ( xi

)

g( xi ) .

i 0

Тогда систему (5.13) можно записать в виде:

c0 (

0 ,

0 )

c1(

0 ,

1 )

...

cm (

0 ,

m )

(

0 , y )

c0 (

1 ,

0 )

c1(

1 ,

1 )

...

cm (

1 ,

m )

(

1 , y )

.

(5.13а)

. . . . . . . . . . . . . . .

c0 (

m ,

0 )

c1(

m ,

1 ) ...

cm ( m , m )

( m , y )

Системы (5.13) и (5.13а) будем называть нормальными системами

уравнений.

Решив

эти

системы,

мы

найдем

коэффициенты c0 ,..., cm , и

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

Алгоритм выбора степени ‘’m’’. В случае, когда m=n мы получим интерполяционный многочлен, поэтому m<<n. Так же необходимо задать

числа	1 и 2, учитывая следующее:
1)	1 >0 и 2>0 должны быть такими, чтобы m находилось между ними;

2)первоначально m выбирают произвольно, но учитывая условие, что m<<n;

3)выбрав m, строят системы (5.13) и (5.13a), решив которые находят

c0 ,..., cm ;

4) используя найденные коэффициенты вычисляется m и проверяется,

попала ли она в промежуток между 1 и 2. Если попала, то степень многочлена выбрана правильно, иначе

57

а) если	m >	1, то степень необходимо уменьшить хотя бы			на
единицу;
б) если	m < 2, то степень необходимо увеличить хотя бы на единицу.
5) затем строить приближающую функцию.
Очень часто для приближения по методу наименьших квадратов
используются	алгебраические		многочлены степени m n, т.е.	k ( x )	xk .
Тогда нормальная система (5.13) принимает следующий вид:
	m	n	n
	(	xij k )c j	yi xik (k= 0,1,…,m).	(5.14)
	j 0 i	o	i 0

Запишем систему (5.14) в развернутом виде в двух наиболее простых случаях m=1 и m=2. В случае многочлена первой степени P1(x)=c0+c1x, нормальная система имеет вид

				n	n
( n 1)c0			(	xi )c1	yi
			i	0	i 0
	n			n	n	(5.15)
(	x )c	0	(	x2 )c	y	x .
	i			i 1		i i
i	0			i 0	i 0

Для многочлена второй степени P2(x)=c0+c1x+c2x2, нормальная система имеет вид

n

n

x2

n

( n 1)c

0

(

x

i

)c (

)c

2

y

i

1

i

i

0

i 0

i 0

n

n

n

n

(

x

i

)c

0

(

x2

)c

(

x3

)c

2

y

i

x

i

.

i

1

i

(5.16)

i

0

i

0

i

0

i

0

n

2

n

x3

n

4

n

x2

(

x

)c

0

(

)c

(

x

)c

2

y

i

i

i

1

i

i

i

0

i

0

i

0

i

0

58

6 Численные методы решения задачи Коши для обыкновенных дифференциальных уравнений и систем дифференциальных уравнений

Будем рассматривать задачу Коши для системы обыкновенных дифференциальных уравнений(ОДУ).Запишем систему в векторной форме

du

f (t, u) ,

(6.1)

dt

где:

u-искомая вектор-функция; t-независимая переменная;

u(t)

(u1 (t),..., um (t)) ; f (t)

( f 1 ,..., f m) , m-порядок

системы;

u (t),...,u

m

(t)

координаты;

t 0; u(0)

u 0 .

1

Запишем систему (6.1) в развернутом виде

d ui

f

(t, u ,..., u

m

) ,

(6.2)

i

dt

1

где:

i=1,...,m;

ui (0) ui0 .

Вслучае i=1 -это будет ОДУ 1-го порядка, а при i=2 - система из двух уравнений первого порядка.

Вслучае i=1 решение задачи Коши предполагает нахождение

интегральной кривой, проходящей через заданную точку и удовлетворяющую заданному начальному условию.

Задача состоит в том, чтобы найти искомую вектор-функцию u, удовлетворяющую (6.1) и заданным начальным условиям.

Известны условия, гарантирующие существование и единственность

решения (6.1) или (6.2).
Предположим, что функции		fi	( i 1,..., m )		непрерывны					по	всем
аргументам в некоторой замкнутой области D={t					a, ui		ui0	b },		где a,b-
известные константы.
Из непрерывности функций следует их ограниченность, т.е. функции
fi сверху ограничены некоторой константой М: | fi					|<M (где М				0) всюду в
области D и пусть в области D функции fi				удовлетворяют условию
Липшица по аргументам u1,..., um . Это значит, что
| f i (t,u 1 ,....,u m )	f i (t,u 1 ,...,u	m )|	L(|u 1	u	1|	....	|u m		u	m|)
для любых двух точек (t, u 1,...., u m)		и	(t, u 1,..., u m)			из области D. Тогда
существует единственное решение задачи (6.1)
u1 u1(t),....,um	um (t) ,определенное при			t	T	min		a, b /M			(6.3)

59

и принимающее при t=0 заданное начальное значение. Существует два класса методов для решения задачи (6.1):

1)семейство одношаговых методов(Рунге-Кутта);

2)семейство многошаговых(m-шаговых) методов.

Сначала рассмотрим одношаговые методы. Для простоты возьмем одно уравнение

				du	f (t, u) ,	(6.4)
				dt
где: u(0) u0 ; t>0.
По оси t	введем равномерную сетку с шагом					, т.е. рассмотрим
систему точек	ωτ tn	n	t,n 0,1,2,..... . Обозначим			через u(t) точное
решение (6.4) , а через		yn	y( tn ) приближенные значения функций u в
заданной системе точек.

6.1 Семейство одношаговых методов решения задачи Коши

6.1.1 Метод Эйлера (частный случай метода Рунге-Кутта)

Уравнение (6.4) заменяется разностным уравнением

yn 1	yn	f (tn , yn ) , n=0,1,2,…, y0	u0 .

В окончательной форме значения yn 1 можно определить по явной формуле

yn 1 yn τ f( tn , yn ).

(6.5)

Вследствие систематического накопления ошибок метод используется редко или используется только для оценки вида интегральной кривой.

Определение 1. Метод сходится к точному решению в некоторой точке t , если yn u(tn ) 0, при , tn t .

Метод сходится на интервале (0,t], если он сходится в любой точке этого интервала.

60