12_116608_1_52628

Пусть на отрезке a,b в узлах сетки заданы значения некоторой

функции f (x), т.е. a x0 x1 x2... xn b, yi f (xi )(i= 0,1,…, n).

Сплайном, соответствующим этим узлам функции 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) многочлен третьей степени

ci

di

Si (x) ai

bi (x xi

)

(x xi )2

(x xi )3 .

(5.8)

2

6

То есть для x xi 1,xi

нужно построить

такую функциюSi ( x),

где

ai ,bi ,ci ,di

подлежат определению. Для всего отрезка интерполирования a,b ,

таким образом, необходимо определить 4 n неизвестных коэффициента.

ci( x

xi )

di

( x

xi )

2

,

S ( x) bi

2

di( x xi ),

S ( x) ci

Si( x) ai

yi .

Доопределим

a0 f (x0) y0 .

Требование

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

функции S(x)

приводит к условиям Si (xi ) Si 1(xi ), (i=0, 1,…,n-1).

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

ai ai 1 bi 1(xi

xi 1 )

ci 1

(xi xi 1 )2

di 1

(xi

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

2

6

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

xi

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

h2

h3

переписать

в

виде

h b

i

c

i

d

f

f

(i=

1,2,…,n).

Из

2

i

i

i

6

i

i

i 1

h2

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

первой

производной

следует

h c

i

d

b

b

(i=

2

i

i

i

i

i 1

2,3,…,n), а

из

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

второй

производной

h

i

d

i

c

i

c

i 1

(i=

2,3,…,n).

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

51

граничными условиями для сплайн-функции в виде S (a) S (b) 0 (концы гибкой линейки свободны).

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

	ci			ci 1 ,c0				cn
hi di									0,( i 1,2,..., n )
			hi2
hi ci					di	bi		bi 1	,( i 2,3,..., n )	(5.9)
	2
	2					3
hi bi			hi		ci		hi	di	fi fi 1 ,( i 1,2,..., n ).
	2					6

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

h c			2(h h	) c h c	6 (	yi 1 yi		yi yi 1	)
	i	i 1	i i 1	i i 1 i 1		hi 1		hi		(5.10)
		cn 0
c0

(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

i

d

i

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

2

6

h

i

i

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)

f( x) Sh( x) M4 h4 ,

8

где M4	max| f IV ( x)|.
	[ a,b]

52

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

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

			n	kx		kx
	многочлена	Tn( x) a0	ak cos		bk sin		, коэффициенты которого
						l
			K 1	l

находятся из системы 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) чтобы сумма квадратов отклонений приближающей функции от таблично заданной была минимальной.

у отклонения

х0 х1 хn-1 хn х

53

Рисунок 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 и
	приравняем их к нулю. Получим систему вида:
		m		n
				2 ( yi c0 0( xi ) ... cm m( xi		)) 0( xi ) 0
		c
				i 0
		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).

54

Определение.	Скалярным произведением функции	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)
	) c1( 1, 1 ) ... cm( 1, m ) ( 1,y)
c0( 1, 0		.
. . . . . . . . . . . . . . .		(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. Если попала, то степень многочлена выбрана правильно, иначе

а) если m > 1, то степень необходимо уменьшить хотя бы на единицу;

б) если m < 2, то степень необходимо увеличить хотя бы на единицу.

5)затем строить приближающую функцию.

Очень часто для приближения по методу наименьших квадратов

используются алгебраические многочлены степени m n, т.е. k(x) xk . Тогда нормальная система (5.13) принимает следующий вид:

m	n	n
( xij k		)cj yixik	(k= 0,1,…,m).	(5.14)
j 0	i o	i 0

55

Запишем систему (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)
( xi )c0 ( xi2 )c1 yixi.
i 0		i 0	i 0

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

n

n

n

(n 1)c0 ( xi )c1 ( xi2 )c2 yi

i 0

i 0

i 0

n

n

n

n

xi )c0

( xi2 )c1 ( xi3 )c2 yixi .

(

i 0

i 0

i 0

i 0

(5.16)

n

2

n

3

n

4

n

2

x

)c

(

x

)c

(

x

)c

y

x

(

i

0

i

1

i

2

i i

i 0

i 0

i 0

i 0

56

6. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧИ КОШИ ДЛЯ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ И СИСТЕМ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

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

		du	f (t,u),
		dt
где: u-искомая вектор-функция; t-независимая переменная;
u(t) (u1(t),...,um (t));	f (t) ( f1,..., f m), m-порядок системы; u1(t),...,

координаты; t 0; u(0) u0 .

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

ddtui f i(t,u1,...,um),

где: i=1,...,m; ui (0) ui0.

(6.1)

um(t)

(6.2)

Вслучае 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)

57

u1 u1(t),....,um um (t) ,определенное при

t T min a,b/M (6.3)

и принимающее при 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], если он сходится в любой точке этого интервала.

58

Определение	2. Метод имеет			р-й порядок точности, если
существует такое число	р>0, для которого			yn u(tn )		O( p ), при	, где:
- шаг интегрирования; O-малая величина порядка p .

	Так как	un 1	un	u (tn) O( ), то метод Эйлера имеет первый

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

6.1.2 Методы Рунге-Кутта

Метод Рунге-Кутта второго порядка точности

Отличительная особенность методов Рунге-Кутта от метода (6.5) заключается в том, что значение правой части уравнения вычисляется не только в точках сетки, но и также в середине отрезков(промежуточных точках).

Предположим, что приближенное значение yn решения задачи в точке

ttn уже известно. Для нахождения yn 1 поступают так:

1)используют схему Эйлера в таком виде

	yn 1 2	yn		f (tn , yn)	(6.6)
	0,5

иотсюда вычисляют y n 1 2 ;

2)воспользуемся разностным уравнением вида

	yn 1 yn	f (tn	0,5 , yn 1 2 ),	(6.7)

откуда найдем значение yn 1. Далее подставим значение yn 1 2 yn 0,5 в

уравнение (6.7). Тогда

yn 1	yn		f(tn 0,5 τ,yn	0,5 τ fn ),	(6.8)
	τ

где fn f (tn , yn ).

Можно показать, что метод (6.8) имеет второй порядок точности, т.е. yn u(t) O( 2).

Метод (6.8) называется методом прогноза и коррекции в том смысле, что на первом этапе решение как бы предсказывается с точностью O( τ ) , а на

втором этапе - с точностью до O( 2 )(второй порядок точности).

59

Будем рассматривать явные методы. Задаем числовые коэффициенты ai , bij , i=2,...,m; j=1,2,...,(m-1) и =1,2,...,m . Последовательно вычисляем функции

k1 f(tn,yn );

k2 f(tn a2τ,yn b21τ k1 );

k3 f(tn a3τ,yn b31τ k1 b32τ k2 );

……………………………………………..

kn f(tn anτ,yn bm1τ k1 ... bmm 1τ km 1).

		yn 1	yn	m
	Затем из формулы			iki	находим значения yn 1. Здесь
				i 1
	ai ,bij , i-числовые параметры,		которые определяются или выбираются из

соображений точности вычислений.

При m=1 и =1 получается метод Эйлера, при m=2 получаем семейство методов

yn 1 yn ( 1k1 2k2),

(6.9)

где: k1 f (tn, yn ); k2 f (tn a2 , yn b21 k1); y0=u0.

Семейство определяет явные методы Рунге-Кутта. Подставив нужные 1 и 2, получаем окончательную формулу. Точность этих методов совпадает с точностью аппроксимирующего метода и равна O( 2 ).

Невязкой, или погрешностью аппроксимации метода (6.9) называется величина

n

un 1 un

1

f (t

n

,u

n

)

2

f (t

n

a ,u

n

b

f (t

n

,u

n

)),

2

21

полученная заменой в (6.9) приближенного решения точным решением.

При 1+ 2=1 получим первый порядок точности. Если же потребовать дополнительно 2b21 2a2 0,5, то получим методы второго порядка точности вида

yn 1 yn

(1 ) f (tn, yn ) f (tn a , yn a f (tn, yn ))

при a 0,5.

Приведем один из методов Рунге-Кутта третьего порядка точности

60