Вычислительная математика. Лекции

А-ортогональную систему можно получить исходя из любой линейно независимой системы векторов fe1; e2; :::; eng:

Следовательно вектор x получается за конечное число итераций, определяемое числом скалярных произведений (b; sk), отличных от нуля. Метод сопряженных направлений, как и метод сопряженных градиентов следует отнести к числу точных методов решения СЛАУ Ax = b:

Замечание. В методе сопряженных градиентов векторы p1; p2; ::: отличные от нуля образуют А-ортогональную систему.

4.6Стационарный многошаговый градиентный метод.

Вградиентном методе векторы fxkg строятся в виде xk+1 = xk krk

(0 < k < 2n )

Тогда

xk+2 = xk krk k+1rk+1 = xk rk k k+1(rk kArk) =

= xk ( k + k+1)rk + k k+1Ark = xk + Ck1rk + Ck2Ark;

xk+3 = xk + Ck1rk + Ck2Ark + Ck3A2rk; :::;

N

X

xk+N = xk + CkiAi 1rk: i=1

28

и

N

X

k xk+N x k2k E + CkiAi kk xk x k :

i=1

NN

Собственные числа операторного полинома E+ P CiAi равны 1+ P Ci ij, где 1; 2; :::; n- собственные

и мы приходим к задаче:

среди полиномов QN ( ) = 1 + C1 + C2 2 + ::: + CN N найти полином, такой что

Такой полином существует и единственен, и мы его построим в следующем параграфе.

4.7Полиномы Чебышева.

Рассмотрим функции 'n(x) = cos(n arccos x); заданные на отрезке [ 1; 1] n = 0; 1; 2; : : :. Ясно, что j'n(x)j 1: При n = 0 '0(x) = 1, при n = 1 '1(x) = x:

Обозначим arccos .

'n+1(x) = cos(n + ) = cos n cos sin n sin ;

'n 1(x) = cos n cos + sin n sin

и

'n+1(x) + 'n 1(x) = 2 cos n cos = 2x'n(x):

Мы получаем реккурентные формулы для вычисления значений функций 'n(x) :

'n+1(x) = 2x'n(x) 'n 1(x):

'2(x) = 2x x 1 = 2x2 1 - четная функция,

'3(x) = 2x(2x2 1) x = 4x3 3x - нечетная функция и т.д.

Таким образом, значения функций 'n(x) могут быть продолжены и для jxj > 1: Тогда 'n(x)- полином степени n с коэффициентом при xn, равным 2n 1: 'n(x) = 2n 1xn + : : : :

Определение. Полиномы Tn(x) = 2n1 1 'n(x) называются полиномами Чебышева. Ясно, что Tn(x) = 1 xn + an 1xn 1 + ::: + a0, и при x 2 [ 1; 1] значения jTn(x)j 2n1 1 :

Нули полинома Чебышева.

Для jxj 1 получаем cos(n arccos x) = 0, n arccos x = 2k2+1 ;

arccos x = 2k2+1 2 , x0k = cos 2kn+1 2 , k = 0; 1; 2; :::(n 1):

На отрезке [ 1; 1] получаем n различных нулей полинома Чебышева. В точке x = 1 значение Tn(1) = 2n1 1 , а Tn( 1) = ( 1)n 1 2n1 1 : Вне отрезка [ 1; 1] нулей полинома Tn(x) нет.

29

Стационарные точки полинома Чебышева.

[x1; 1] полином Rn 1 имеет не менее одного нуля, на отрезке [x2; x1] полином Rn 1 имеет не менее одного

нуля и на каждом отрезке [xk; xk 1], k = 2; 3; :::; (n 1); а также на отрезке [ 1; xn 1] полином имеет не менее одного нуля. Тогда на всем отрезке [ 1; 1] полином Rn 1(x) имеет не менее n нулей, что возможно

только если Rn 1(x) 0 на (1; 1).

Если же в некоторой точке xkTn(xk) = Qn(xk), то в этой точке Q0n(xk) = Tn0 (xk) и xk есть ноль, порядок которого не меньше двух. С учетом кратности полином Rn 1 имеет на отрезке [ 1; 1] не менее n нулей и

поэтому Tn(x) = Qn(x):

Полином Чебышева - полином наименее отклоняющийся от нуля на отрезке [ 1; 1] среди всех полиномов вида xn + cn 1xn 1 + ::: + c0:

Построим полином n( ) степени n наименее отклоняющийся от нуля на отрезке [a; b] среди всех полиномов вида 1 + c1 + c2 2 + ::: + cn n; предполагая что a > 0: Линейное преобразование

отображает отрезок a b на отрезок [ 1; 1]: Полином наименее отклоняющийся от нуля на отрезке

Заметим, что bb+aa < 1 и значение Tn( bb+aa ) следует вычислять, исходя из полиномиального представления полинома Чебышева.

4.8Многошаговый стационарный градиентный метод

(продолжение).

Поставленная в конце x6 задача имеет единственное решение:

При реализации метода достаточно знать только нули 0k этого полинома:

30

Обозначим

1

k = k ; N ( ) = (1 0)(1 1) : : : (1 N 1)

и

N (A) = (E 0A)(E 1A) : : : (E N 1A):

Построим последовательность векторов xk+i :

xk+1 = xk 0rk и xk+1 x = xk x 0A(xk x ) = (E 0A)(xk x ); xk+2 = xk+1 1rk+1 и xk+2 x = (E 1A)(E 0A)(xk x )

: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :

xk+N = xk+N 1 N 1rk+N 1 и xk+N x = (E N 1A)(E N 2A) : : : (E 1A)(E 0A)(xk x ) =N (A)(xk x ):

Ясно, что порядок перечисления шагов при построении xk+N безразличен. Для k N (A) k2 верна оценка

pp !N

Найдя xN+k повторяет цикл, возможно с другим набором шагов , изменяя степень полинома Чебышева.

Заключение. Кроме методов решения СЛАУ и минимизации значений квадратичной функции существуют и другие методы (метод вращений, метод покоординатного спуска и др.) Для минимизации значений произвольных (дифференцируемых) функций многих переменных применяется метод последовательной линеаризации в окрестности точек xk с последующей минимизацией значений получаемых квадратичных функций.

31

Глава 5

ИНТЕРПОЛИРОВАНИЕ ФУНКЦИЙ.

ВВЕДЕНИЕ. Общая постановка задачи.

Рассматривается линейное метрическое пространство F . В этом пространстве выбрано (n+1) линейно независимых ("базисных") элементов '0; '1; :::; 'n и задано (m+1) функционалов J0; J1; :::; Jm.

Задача. Для любого элемента f 2 F требуется построить линейную комбинацию f~ элементов 'n

n

f~ = XCk'k;

k=0

такую что Ji(f) = Ji(f~); i = 0; 1; :::; m.

В общем случае элементы f и f~ не совпадают. Но мы не можем их различить, зная только (m+1) значений функционалов J0; J1; :::; Jm на этих элементах. Как правило построение элемента f~ не является конечной целью исследований. Простая структура элемента f~ обеспечивает простоту выполнения дальнейших операций. Поэтому для оценки точных результатов необходимо уметь оценивать близость элементов f и f~, т.е. оценивать расстояние (f; f~) в метрике пространства F .

Поставленная задача является классической задачей вычислительной математики, она изучалась многими учеными, начиная с И.Ньютона (1711г.) и Ж.Лагранжа (1806г.). Существует множество учебников и учебных пособий по этой тематике, среди них отметим:

-Н.Бахвалов, Н.Жидков, Г.Кобельков. Численные методы. 2003г. -В.Вержбицкий. Основы численных методов. 2002г.

В этих учебниках приведена обширная библиография.

И н т е р п о л и р о в а н и е ф у н к ц и й о д н о й п е р е м е н н о й.

5.1Интерполирование функций алгебраическими полиномами.

В качестве пространства F задается пространство [a; b] функций заданных и непрерывных на отрезке

[a; b]:

(f; y) = max jf(x) f(y)j

x2[a;b]

"Базисными"элементами выбраны функции

'k(x) = xk; k = 0; 1; :::; n; x [a; b]

На отрезке [a,b] заданы различные точки xi (узлы), а в качестве функционалов Ji(f) рассматриваются функционалы

Ji(f) = f(xi); i = 0; 1; :::; n; (n = m)

Линейные комбинации "базисных"элементов-полиномы степени n, рассматриваемые на [a,b].

Задача. Для заданной непрерывной на [a,b] функции f(x) построить полином f~n(x) степени n такой что

n

X

f~n(x) = Ck(f)xki = f(xi)

k=0

32

во всех узлах xi

Решение этой задачи легко получить: достаточно решить систему линейных алгебраических уравнений относительно неизвестных значений C0; C1; :::; Cn:

n

X

Ck(f)xki = fi; i = 0; 1; :::; n; fi = f(xi):

k=0

Матрица этой системы

и её определитель - определитель Вандермонда-отличен от нуля. По формулам Крамера

где Aik-алгебраические дополнения элементов k-того столбца определителя . Тогда

i=0

Полиномы n-ой степени li(x) зависят только от выбора узлов x0; :::; xn. Полином f~n(x) называется интерполяционным полиномом для функции f(x) с узлами x0; :::; xn

5.2Интерполяционный полином в форме Лагранжа.

Значение полиномов li содержит вычисление определителей порядка (n + 1). Так как эти определители

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

n

полиномы li без непосредственного вычисления определителей. Равенства f(xj) = f~n(xj) = P f(xi)li(xj)

i=0

будут выполнены, если потребовать, чтобы

Единственным полиномом степени n, обладающим этим свойством, является полином

(x x0)(x x1):::(x xi 1)(x xi+1):::(x xn) (xi x0)(xi x1):::(xi xi 1)(xi xi+1):::(xi xn)

Компактную запись полинома li(x) можно получить, если ввести обозначение !n+1(x) = (xi x0)(xi

x1):::(x(i) xi 1)(xi xi+1):::(xi xn). Ясно, что производная функции !n+1 в точке x = xi равна !n0 +1(xi) = (xi x0)(xi x1):::(x(i) xi 1)(xi xi+1):::(xi xn), и тогда li(x) = !n+1(x) Полученная таким образом

форма интерполяционного полинома обозначается Ln(f; x):

n

X

f~n(x) = Ln(f; x) = i=0 f(xi)(x xi)!n0 +1(x)

33

и называется интерполяционным полиномом в форме Лагранжа. Отметим три свойства интерполяционных полиномов, не зависящих от конкретной формы интерполяционного полинома:

1. Интерполяционный полином не зависит от порядка перечисления узлов.

достаточно заметить, что коэффициенты полиномов li(x) при xn равны 1. Особенностью представления интерполяционного полинома в форме Лагранжа является возможность быстрого построения интерполяционного полинома для любой функции f F при фиксированном наборе узлов, если полиномы li(x) уже найдены. Однако при добавлении нового узла придется заново строить полиномы li(x).

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

и эта величина не зависит от порядка перечисления узлов. Проведя элементарные выкладки, легко показать (по индукции), что разделенные разности k-го порядка равны

k

X f(xi) f(x0; x1; :::; xk) = i=0 !k0 +1(xi)

и их значения не зависят от порядка перечисления узлов в наборе x0; x1:::; xk: Согласно свойству 3) коэффициент ak при xk в интерполяционном полиноме f~k(x) равен

ak = f(x0; x1; :::; xk):

34

5.4Интерполяционный полином в форме Ньютона.

Пусть f~k(x)-интерполяционный полином для функции f(x), построенный по узлам x0; x1; :::; xk 1. Тогда разность f~k(x) f~k 1(x)-полином степени k, обращающийся в ноль в точках x0; x1; :::; xk 1. Ясно, что

f~k(x) f~k 1(x) = a(x x0)(x x1):::(x xk 1) = a!k(x);

где a равен коэффициенту полинома f~k(x) при xk:

a = ak = f(x0; x1; :::; xk):

Интерполяционнай полином для функции f(x), построенный по одному узлу x0 равен

f~0(x) = f(x0)

Поэтому естественно определить разделенную разность нулевого порядка как f(x0), а значение !0(x) 1. Ясно что

f~n(x) = f~0(x) + [f~1(x) f~0(x)] + [f~2(x) f~1(x)] + :::+

Такая форма представления интерполяционного полинома f~n(x) обозначается Nn(f; x) и говорят, что интерполяционный полином представлен в форме Ньютона:

n

X

f~n(x) = Ln(f; x) = Nn(f; x) = f(x0; x1; :::; xk)!k(x):

k=0

При введении нового узла xn+1

Nn+1(f; x) = Nn(f; x) + f(x0; x1; :::; xn; xn+1)!n+1(x):

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

f(xn; xn+1); f(xn 1; xn; xn+1); :::; f(x1; x2; :::; xn; xn+1); f(x0; x1; :::; xn; xn+1):

5.5Оценка методической погрешности.

Обозначим rn(f; x) = f(x) f~n(x) и оценим

max jrn(f; x)j = kf f~nkC[a;b]

x2[a;b]

для заданной функции f(x) и для фиксированного набора узлов x0; x1; :::; xn. Для краткости будем обозначать rn(x) = rn(f; x). Ясно, что без дополнительных ограничений на функцию f(x) невозможно получить оценку величины rn(x). Будем предполагать, что функция f(x) имеет на [a,b] непрерывную производную

порядка (n+1) и известна оценка:

jfn+1(x)j Mn+1:

Пусть x -произвольная точка отрезка [a,b]. Оценим rn(x ). Ясно, что если x совпадает с одним из узлов xi, то rn(x ) = 0. Поэтому будем рассматривать точку x не совпадающую ни с одним из узлов. Построим вспомогательную функцию (x):

(x) = rn(x) !n+1(x) rn(x )

!n+1(x )

Эта функция имеет непрерывные производные до порядка (n+1) включительно. В точках x0; x1; :::; xn и в точке x функция (x) обращается в ноль. Таким образом, эта функция имеет на [a,b] не меньше (n+2)

35

нулей. Тогда производная 0(x) имеет не менее (n+1) нулей и n+1(x) имеет на [a,b] по крайней мере один ноль. Следовательно существует точка [a; b], в которой n+1( ) = 0. Так как

Эту оценку нельзя улучшить:взяв функцию f(x):

для которой f(n+1)(x) = Mn+1, мы получим

1

jrn(f; x)j = (n + 1)!Mn+1j!n+1(x)j:

5.6Выбор узлов интерполирования.

Рассмотрим постановку задачи интерполирования, в которой узлы x0; x1; : : : ; xn не фиксированы и их можно выбирать для улучшения оценки методической погрешности. Исходя из оценки (2) будем выбирать узлы xi из условия

min ( max j!n+1(x)j):

fx0;x1;:::;xng x2[a;b]

Так как !n+1(x) - полином степени (n+1) и !n+1(x) = (x x0)(x x1):::(x xn) = 1 xn+1+anxn+:::+a1x+a0, то оптимальный полином - полином, наименее отклоняющийся от нуля на отрезке [a; b]. Если [a; b] = [ 1; 1], то !n+1(x)-полином Чебышева Tn+1(x), нулями которого являются числа xi, равные

2i + 1

xi = cos 2n + 1 ; i = 0; 1; :::; n:

Для построения полинома !n+1(x) выберем узлы xi - нули полинома наименее отклоняющегося от нуля на отрезке [a; b]:

2i + 1

xi = cos 2n + 2 ; i = 0; 1; :::; n:

Согласно оценке значений полинома, наименее отклоняющегося от нуля, получаем

36

5.7Разделенные разности и производные функции.

Кнабору узлов x0; x1; :::; xn добавим еще один узел x 2 [a; b]. Тогда

Nn+1(f; x) = Nn(f; x) + f(x0; x1; :::; xn; x )!n+1(x)

и в точке x = x

Nn+1(f; x ) = f(x ) = Nn(f; x ) + f(x0; x1; :::; xn; x )!n+1(x ):

Так как x - произвольная точка отличная от xi, то получаем тождество

rn(f; x) = f(x0; x1; :::; xn; x)!n+1(x);

верное и для всех значений x 2 [a; b]. Согласно равенству (1)

где 2 [a; b]:

5.8Численное дифференцирование. Методическая погрешность

численного дифференцирования.

Построив интерполяционный полином, легко находить приближенные значения интегралов и производных исходной функции f(x). В этом параграфе мы оценим точность приближения производных f(m)(x) при замене функции f(x) интерполяционным полиномом в предположении, что функция f имеет [a; b] непрерывные производные до порядка n + m + 1 включительно. При этом мы будем представлять интерполяционный полином в форме Ньютона. Ясно что m должно быть не больше n : m n.

Представим rn(x) в виде rn(x) = f(x0; x1; :::; xn; x)!n+1(x). Ясно что f(x0; x1; :::; xn; x) имеет столько же непрерывных производных сколько и функция f(x).

X

rn;m(f; x) = Ck dx f(x0; :::; xn; x)dx(m k) !n+1(x):

m k

k=0

Мы получим оценку rn;m(f; x), если сумеем оценить производные разделенных разностей

dk

dxk f(x0; :::; xn; x):

Для этого рассмотрим произвольную функцию g(x) 2 Ck[a; b] и её разделенную разность g(x; x + h; x + 2h; :::; x + kh) порядка k. Согласно (3) её можно представить в виде

g(x; x + h; x + 2h; :::; x + kh) = k1!g(k)(x);

37