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

и

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

h!0

Рассмотрим в качестве g(x) функцию g(x) = f(x0; x1; :::; xn; x) - разделенную разность порядка (n + 1) функции f(x).

тельность h ! 0 и последовательность h. Последовательность h можно и не иметь предела, но так как она ограничена, то из неё можно выбрать последовательность, сходящуюся к некоторой точке [a; b]. Для этой подпоследовательности существует предел

Замечание: Для равностоящих узлов xi; x0 = a; xn = b можно получить оценку rn;m(f; x) для функций f Cn+m[a; b].

r1;1(f; x) = f(x0; x1; x)dxd !2(x) + f0(x0; x1; x)!2(x) = 2!1 f00( 1)(2x h) + 3!1 f000( 2)(x x0)(x x1)

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

r1;1(f; x) = f0(x) [f(x0) + f(x1) f(x0)(x x0)]0 = f0(x) f0( ) h

и

jr1;1(f; x)j M2h:

38

5.9Полная погрешность численного дифференцирования.

Как правило значение функции f или вследствие погрешности их вычисления или вследствие их измерений известны с некоторой точностью, т.е. известны значения f (x) = f(x) + (x), где (x) - погрешность задания функции f(x).

Тогда полная погрешность интерполирования равна

f(x) (Ln(f; x) + Ln( ; x)) = rn(f; x) Ln( ; x):

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

В этом параграфе мы оценим полную погрешность численного дифференцирования. Предполагая что функция f 2 Cn+2[a; b], получим

Относительно функции (x) мы можем предполагать только её ограниченность j (x)j :

Остановимся на одном частном случае.

Рассмотрим два узла x0 и x1 = x0 +h и функцию f 2 C2[x0; x1]. Методическая погрешность численного дифференцирования оценивается jr1;1(f; x)j M2h, а так как

Ln( ; x) = (x0) + (x1) (x0)(x x0) h

мальным. Ясно, что минимум достигается при

r 2 p h = :

M2

p p

При таком выборе h полная погрешность численного дифференцирования не превосходят 2 2M2 . Уменьшая, например, погрешность задания функции f(x) в 100 раз мы получим уменьшение полной погрешности численного дифференцирования всего в 10 раз при одновременном уменьшении шага в 10 раз.

5.10Интерполирование с равноотстоящими узлами.

Рассматриваются узлы a = x0; xi = x0 + ih; xn = x0 + nh = b и значения fi = f(xi) непрерывной на [a; b] функции f(x). Определим конечные разности kfi k-того порядка:

1f0 = f1 f0

39

k+1f0 = kf1 kf0

и их связь с разделенными разностями f(x0; x1; :::; xk) порядка k.

Ясно, что f(x0; x1) = 1fo . Легко показать, что верна формула

h

Cn+1xn+1 + + C1x + C0

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

Таким образом для оценки rn(f; x) следует оценить j!n+1(x)j: Из последнего представления интерполяционного полинома

n

f~n(x) = X kf0 !k(x) k!hk

k=0

можно получить интересное свойство f~n(x). Пусть в процессе вычислений получены значения kf0; k = 0; 1; : : : ; n. Оценим вклад k-того слагаемого

При малых t величина t(1 t)k1 . Таким образом, вычислив значения 4kf0 можно судить о влиянии k-того слагаемого в интерполяционной формуле на значения интерполяционного полинома для x 2 [a; a + h]. Это влияние убывает с ростом k достаточно медленно. Немецкий математик Рунге около 1900 г. рассмотрел на отрезке [ 1; 1] бесконечно дифференцируемую функцию

f(x) =

1

1 + 25x2

40

и показал, что интерполяционный полином f~n(x), построенный по равноотстоящим узлам, не стремится с ростом n к этой функции:

при n ! 1: Завершая этот параграф, отметим, что можно строить интерполяционный полином, взяв в качестве начальной точки значение x = b и шаг h. Можно взять в качестве x0 середину отрезка [a; b] и одновременно брать узлы x0 + kh; x0 kh. Эти и другие возможности приведены в книге В.М.Вержбицкого "Основы численных методов".

5.11Интерполирование с кратными узлами (интерполяционный

полином Эрмита).

Далее мы рассмотрим другие постановки задач интерполирования. Пусть в узлах x0; x1:::xn помимо значений f(xi) известны (вычислены,измерены) значения производных функции f(x):

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

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

Обозначим N = maxfmfigg: Для каждого узла xi рассматриваются все производные до порядка mi включительно. Число всех величин в узлах равно

m0 + m1 + ::: + mn + n + 1 = m + 1; m = n + N:

Требуется построить полином, удовлетворяющий всем поставленным условиям. Такой полином должен иметь степень m и этот полином, если он существует, определяется единственным образом. Действительно, если разность двух таким полиномов Pm(x) Qm(x) = Rm(x) не равна тождественно нулю, то степень полинома Rm не превосходит m. Но узлы xi являются его нулями кратности mi:

Rm(x) = a(x x0)m0+1(x x1)m1+1:::(x xn)mn+1:

Стало быть степень полинома Rm(x) равна m + 1, что невозможно.

Для искомого полинома можно построить систему линейных алгебраических уравнений относительно его неизвестных коэффициентов. Она имеет единственное решение. Таким образом поставленная задача имеет решение и это решение единственно.

Полином, удовлетворяющий условиям задачи, называется интерполяционным полиномом Эрмита и обозначается Hm(f; x).

Укажем "наивный"способ построения полинома Hm(x). В непересекающихся h- окрестностях каждого узла xi выберем mi дополнительных точек xi;1; xi;2; :::; xi;mi . Вычислим значения в этих точках, используя формулу Тейлора для функции f(x) с известными значениями f(xi); f0(xi); :::; f(mi)(xi).

Вместе с узлами x0; x1; :::; xn мы получим ещё m0+m1+:::+mn узлов, в которых значения функции f(x) известны приближенно. Построим интерполяционный полином Lm(f; h; x) по этим узлам. Тогда, устремив h к нулю, получим

Hm(f; x) = lim Lm(f; h; x):

h!0

Ясно, что такой метод неконструктивен. Но он позволяет получить представление rm(f; x) = f(x)

41

5.12Метод построения интерполяционного полинома Эрмита.

Укажем один из возможных методов построения полинома Hm(f; x).

Функция F1(x) неизвестна. Однако её значения и значения всех производных, участвующих в таблице, можно вычислить. Действительно

dxd [f(x) L(f; x)] = f0(x) L0(f; x) = F1(x)!0(1; x) + F10!(1; x);

ив узлах таблицы f0(xi) L0(f; xi) = F1(xi)!0(1; xi); F1(xi) = f0(xi)0 L0(f;xi) : Далее

!(1;xi)

d22 [f(x) L(f; x)] = f00(x) L00(f; x) = F1(x)!00(1; x) + 2F10!0(1; x) + F100!(1; x): dx

В тех узлах, где заданы значения вторых производных, получаем значения F10(xi):

F10(xi) = f00(xi) L00(f; xi) F1(xi)!00(1; xi)

2!00(1; xi)

Продолжая дифференцирование функции f(x) L(f; x), получаем значения функции F1(x) и всех её производных в узлах, входящих в таблицу.

Таким образом для функции F1(x) получаем задачу кратного интерполирования. Ясно что

f(x) = L(f; x) + !(1; x)F1(x)

3. Построим интерполяционный полином L(F1; x) по значениям F1(x) в узлах и представим разность

F1(x) L(F1; x) в виде

F1(x) L(F1; x) = !(2; x)F2(x)

Тогда

f(x) = L(f; x) + !(1; x)L(F1; x) + !(1; x)!(2; x)F2(x):

Аналогично пункту (2) получим задачу кратного интерполирования для функции F2(x):

4. Продолжая этот процесс, получим задачу кратного интерполирования для функции FN 1(x), построим интерполяционный полином L(FN 1; x) и разность FN 1(x) L(FN 1; x) представим в виде

FN 1(x) L(FN 1; x) = !(N; x)FN (x):

Для функции FN (x) получим только её значение в узлах исходной таблицы, в которых заданы значения производных f(x) наибольшего порядка N. По значениям FN (x) в этих узлах построим интерполяционный полином L(FN ; x) и FN (x) = L(FN ; x) + r(FN ; x) Окончательно получаем

f(x) = L(f; x) + !(1; x)L(F1; x) + :: + !(1; x)!(2; x)::!(N; x)[L(FN ; x) + r(FN ; x)]

Пренебрегая величиной r(FN ; x), получим интерполяционный полином Эрмита:

Hm(f; x) = L(f; x) + !(1; x)L(F1; x) + ::: + !(1; x)!(2; x):::!(N; x)L(FN ; x)

Представление f(x) Hm(f; x) можно получить, минуя представление r(FN ; x). Согласно (5):

Замечание: Если значение функции f и всех её производных до порядка n включительно заданы только в одном узле, то

Hn(f; x) = f(x0) + f0(x0)(x x0) + 2!1 f00(x0)(x x0)2 + ::: + n1!fn(x0)(x x0)n

42

n = 2; m = 5; N = max mi = 2:

Построить интерполяционный полином Эрмита H5(f; x), удовлетворяющий этой таблице значений.

1.По значениям f(x) в узлах -1,0,1 построим интерполяционный полином L(f; x) : L(f; x) = 1 + x.

2.Разность f(x) L(f; x) представим в виде f(x) L(f; x) = !(1; x)F1(x), где (1; x) = !(x x0)(x x1)(x x2) = x3 x, F1(x)-неизвестная функция.

3. Для функции F1(x) построим интерполяционный полином по её значениям в узлах x0; x1 : L(F1; x) = 1 x и представим разность F1(x) L(F1; x) в виде F1(x) L(F1; x) = !(2; x)F2(x) и

Производная функции F1(x) равна

F10(x) = 1 + (2x + 1)F2(x) + !(2; x)F20(x)

4. Для функции F2(x) строим интерполяционный полином по её значению с точке x0 : L(F2; x) = 1 и

F2(x) = 1 + r(F2; x):

Если F2(x) дифференцируемая функция, то r(F2; x) = F20( )(x x0). Итак

f(x) = 1 + x + !(1; x)(1 x + (x + 1)x(1 + r(F2; x))) = 1 + x + (x3 x)(1 x + x2 + x + x(x + 1)r(F2; x)) = = 1 + x5 + (x 1)(x + 1)2x2r(F2; x):

Пренебрегая величиной r(F2; x), получаем

H5(f; x) = 1 + x5:

Согласно (5)

f(x) H5(f; x) = (x + 1)3x2(x 1)6!1 f(6)( );

если f 2 C6[x0; x2].

43

5.13Интерполяционный процесс.

Рассмотрим последовательность разбиений n отрезка [a,b] узлами x(1n); x(2n); :::; x(nn), требуя чтобы с ростом n диаметр разбиения стремился к нулю.

Для функции f(x) построим последовательность интерполяционных полиномов Ln( n; f; x). Разбиениеn задает интерполяционный процесс ( n; f) для функции f. Обозначим r( n; f; x) = f(x) Ln( n; f; x).

Говорят, что интерполяционный процесс ( n; f) сходится в точке x 2 [a; b], если lim r( n; f; x ) = 0 при n ! 1.

Говорят, что интерполяционный процесс сходится равномерно на [a; b], если max jr( n; f; x)j ! 0 при

x2[a;b]

n ! 1. В этом случае пишут r( n; f; x) 0 при n ! 1. Теорема Марцинкевича.

Для каждой непрерывной функции f(x) существует такой интерполяционный процесс, что r( n; f; x) 0. В этой теореме n = n(f).

Теорема Фабера.

Для любой последовательности разбиений n отрезка [a; b] существует непрерывная функция f(x),

такая что max jr( n; f; x)j не стремится к нулю при n ! 1.

x2[a;b]

В этой теореме f = f( n).

Для отрезка [ 1; 1] при разбиении разноотстоящими узлами пример такой функции привел Рунге:

1

f(x) = 1 + 25x2 :

Для функции f(x) = jxj для такого же разбиения отрезка [ 1; 1] Бернштейн показал, что lim r( n; f; x) =

h!1

0 только в точках x = 1; x = 0; x = 1.

Рассмотренные постановки задач интерполирования имеют существенные недостатки: -увеличение степени полиномов при увеличении числа узлов.

-величина rn(f; x) = f(x) f~n(x) оценивается через значения производной функции f порядка (n + 1) в предположении, что f 2 Cn+1[a; b].

-увеличение числа узлов не гарантирует в общем случае уменьшения max jrn(f; x)j.

x2[a;b]

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

5.14Сплайны класса Sn;k. (Кусочно-полиномиальное

интерполирование)

Пусть задана возрастающая последовательность узлов

x0 < x1 < x2 < ::: < xN 1 < xN

Функция S(x) 2 Ck[x0; xN ] называется сплайном класса Sn;k, если на каждом отрезке между узлами она

является полиномом степени не выше n.

Если k = n, то Sn;n-множество степени не выше n, заданных на [x0:::xN ]. Величина =. n k 1 называется дефектом сплайна.

Для фиксированной последовательности узлов xi множество функций класса Sn;k является линейным подпространством пространства Ck[x0; xN ]. Покажем, что это подпространство конечномерно, построив базис в классе Sn;k.

Обозначим

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

44

Ясно, что эти базисные функции 'i;j(x) линейно независимы. Назовем их исходными базисными функциями класса Sn;k.

В каждом внутреннем узле xi значения 'i;j(x); '0i;j(x); :::; 'ki;j(x) равны нулю, поэтому линейные комбинации этих базисных функций являются сплайнами класса Sn;k. С другой стороны каждый сплайн этого класса можно представить в виде линейной комбинации исходных базисных функций.

Введем сквозную нумерацию исходных базисных функций:

Постановка задачи интерполяции. Для функции f 2 C[x0; xN ] построить сплайн S(x) класса Sn;k, такой

что S(xi) = f(xi); i = 0; 1; :::; N:

@

P

Представляя искомый сплайн в виде S(x) = Cj; 'j, приходим к задаче определения коэффициентов

j=1

Cj.

Таким образом задача интерполяции кусочно-полиномиальными функциями является частным случаем задачи интерполирования. В общем случае число условий N + 1 меньше @ и коэффициенты Cj определяются неоднозначно. Решение поставленной задачи можно упростить, если рассмотреть новый базис si(x); i = 0; 1; :::; N, где каждая функция si(x) является линейной комбинацией исходных базисных функций.

Например, для класса S1;0 такими линейными комбинациями являются сплайны класса S1;0, построенные по узлам xi 1; xi; xi+1, равные нулю при x xi 1 и при x xi+1, а в узле xi значение si(xi) = 1.

5.15Построение интерполяционного сплайна класса S3;2.

По сравнению со сплайнами классов S3;0 и S3;1 сплайны класса S3;2 имеют наибольшую гладкость. Необходимость построения интерполяционных сплайнов класса S3;2 возникает при решении многих

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

ZxN d4S 2

dx4

dx

x0

Естественно возникает задача минимизации P (S) при условиях S(xi) = f(xi) и условиях непрерывности S(x); S0(x); S00(x). Эта задача-классическая задача вариационного исчисления. Оптимальная функция S(x) должна удовлетворять уравнению Эйлера

ZxN d4S 2

dx4

dx = 0

x0

и граничным условиям S00(x0) = 0; S00(xN ) = 0. Следовательно S(x) является полиномом степени 3 на каждом отрезке [xi 1; xi], а на всем отрезке [x0; xN ]-сплайн класса S3;2 с дополнительными граничными условиями. Обозначим Si(x) значения сплайна S(x) при x 2 [xi; xi+1], а через hi обозначим длину отрезка

[xi; xi+1]; hi = xi+1 xi. Так как Si(xi) = f(xi), то

Si(x) = fi + bi(x xi) + 12Ci(x xi)2 + 3!1 di(x xi)3;

где коэффициенты bi; Ci; di подлежат определению.

1. Условие непрерывности S(x) во внутренних точках xi:

45

3. Условие непрерывности S00(x) во внутренних узлах xi:

Si00(xi) = Si00 1(xi); Ci = Ci 1 + di 1hi 1

или

di 1 = Ci Ci 1 ; i = 1; 2; :::; N 1::::::::::(C)

Таким образом для 3N неизвестных коэффициентов мы получили N + 2(N 1) = 3N 2 уравнений. В теории упругости для балки с фиксированными торцами условия равновесия имеют вид

S00(x0) S0(x0) = 0

S00(xN ) S0(xN ) = 0;

где характеризует упругие свойства внешней среды. Эти два соотношения дают необходимые условия для разрешимости системы линейных уравнений.

В случае, если к торцам балки не приложены силы, то коэффициент = 0 и для граничных точек получаем: S00(x0), т.е. C0 = 0 и S00(xN ) = 0, т.е. SN00 1(xN ) = 0; CN 1 + dN 1hN 1 = 0: Таким образом получаем два недостающих условия:

C0 = 0; dN 1 = CN 1 hN 1

Последнее условие можно включить в условия (С) при i = N, положив формально CN = 0 : dN 1 =

CN CN 1 .

hN 1

Метод решения полученной системы состоит в последовательном исключении неизвестных di, затем bi.

Врезультате мы получим систему линейных уравнений для определения коэффициентов Ci. Условия (С):

46

Умножая эти уравнения на 6, получаем систему линейных алгебраических уравнений с матрицей размерности (N 1) (N 1):

и с вектор-столбцом правой части:

Эта матрица трёхдиагональная с диагональным преобладанием. Решение системы уравнений можно получить методом прогонки, прогоночные коэффициенты по модулю меньше 1, метод устойчив. Для других граничных условий изменяются первое и последнее уравнения, и присоединяются еще два уравнения при i = 0 и при i = N.

5.16 Оценка методической и полной погрешности интерполяционного чертежного сплайна.

Для этих оценок рассмотрим равномерную сетку узлов: hi = h. Система уравнений для коэффициентов

47