Schisla_1

	0.8
	0.6
	0.4
	0.2
	0	0.2	0.4	x	0.6	0.8	1
	–0.2
	–0.4
	–0.6
	–0.8

Рис. 3: Пример сильно осциллирующей функции.

Выгоднее (если возможно) представить подынтегральную функцию в виде

f(x, ω) = Real(φ(x, ω)ei!x),

или иначе

f(x, ω) = φ1(x, ω) sin ωx + φ2(x, ω)cosωx,

где функции φ′1(x, ω) и φ′2(x, ω) имеют значительно меньше перемен знака на [a, b].

Тогда для I(f) получится выражение

∫b ∫b ∫b

I(f) = f(x, ω)dx = ρ1(x)φ1(x, ω)dx + ρ2(x)φ2(x, ω)dx,

a a a

где ρ1(x) = sin ωx и ρ2(x) = cos ωx можно рассматривать как весовые функции.

Пример. Вывод квадратурной формулы с весом ρ(x) = sin ωx и двумя узлами x = a и x = b.

Для вывода формулы вычисления приближенного значения интеграла

∫b

I(f) = ρ(x)f(x, ω)dx.

a

61

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

f(x, ω) = L2(x, ω) + r2(x, ω),

где L2(x, ω) = f(a, ω)

x − b

+ f(b, ω)

x − a

– интерполяционный многочлен

a − b

b − a

f′′(ξ, ω)

Лагранжа, r2(x, ω) =

(x − b)(x − a) – его погрешность, ξ (a, b).

2

Теперь

I(f) = ∫a b ρ(x)f(x, ω)dx ≈ S2(f) = ∫a b ρ(x)L2(x, ω)dx =

= C1f(a, ω) + C2f(b, ω),

b

x − b

b

x − a

где

C

=

∫a

sin ωx dx

,

C

=

∫a

sin ωx dx

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

1

a − b

2

b − a

квадратурной формулы. Вычислить их можно, интегрируя по частям,

−

x

b) cos ωx

sin ωx

b

cos ωa

−

1

(−

(

) a=

sin ωb

sin ωa

C1 =

−

+

−

−

,

a b

ω

ω2

ω

ω2(b

a)

аналогично

C

2

=

−

cos ωb + sin ωb − sin ωa .

ω

ω2(b − a)

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

∫b

R2(f)(f) = I2(f) − S2(f) = (f(x, ω) − L2(x, ω))dx =

a

∫b

= sin ωx f′′(ξ, ω)(x − b)(x − a)dx,

2

a

из которого следует оценка

|R(f)| 6 M2(b − a)3/6,

где M2 = max |f′′(x)|.

a6x6b

Если аналогичным образом вывести квадратурную формулу для весовой функции ρ(x) = ei!x, используя интерполяционный многочлен, построенный по трем узлам, то получится формула Филона.

62

4.7 Вопросы и упражнения для самоконтроля.

Вопросы

1.Как оценивается погрешность квадратурной формулы Ньютона - Котеса, построенной по 5 узлам?

2.Что такое алгебраическая точность квадратурной формулы? Какова она у квадратурной формулы Ньютона-Котеса, построенной по четырем узлам?

3.Какова алгебраическая точность квадратурной формулы Симпсона?

4.Какова алгебраическая точность квадратурной формулы Гаусса, построенной по 3 узлам?

5.Как оценивается погрешность квадратурной формулы Гаусса, построенной по 3 узлам?

6.Какие многочлены называются ортогональными с весом ρ(x)?

7.Какой вид имеет рекурентная формула для построения системы ортогональных с весом ρ(x) многочленов?

8.Для каких весовых функций ρ(x) существует система ортогональных многочленов?

9.Какими свойствами обладают корни ортогональных многочленов?

10.Как экономичнее вычислять интегралы от сильно осциллирующих функций?

5Численное дифференцирование

Формулы численного дифференцирования функции, заданной таблич-

но, можно построить также, как и квадратурные.

63

5.1Вывод формул численного дифференцирования с помощью интерполяции

Пусть значения функции f(x) заданы на [a, b] в узлах сетки x0, ..., xn−1. Приблизим функцию интерполяционным многочленом, построенным по этим узлам

f(x) = Ln(x) + rn(x),

где

	n−1				n−1		−
	∑k				∏̸	x
Ln(x) =		f(xk)				xk	− xj		,
	=0				j=0;j=k			xj
а rn(x) = f(x)−Ln(x) =	f(n)(ξ(x)) n−1			(x−xk) – погрешность интерполяци-
	n!		=0
			∏k							(k)	(x)
онного многочлена Лагранжа, ξ(x) (a, b). Тогда для производной f

производной получается

f(k)(x) = L(nk)(x) + rn(k)(x).

Отсюда, положив f(k)(x) ≈ L(nk)(x), получим расчетную формулу для численного дифференцирования, а также оценку ее погрешности

f(k)(x) − L(nk)(x) = rn(k)(x).

Выводить формулы численного дифференцирования и обсуждать погрешность выведенных формул для произвольного числа узлов n громозд-

ко. Обсудим все это на примере n = 3.

Построим формулы численного дифференцирования в случае равноотстоящих узлов xi = x0 + ih, i = 0, 1, 2. Запишем интерполяционный много-

член в виде

L3(x0

+ th) = y0

(t − 1)(t − 2)

+ y1

t(t − 2)

+ y2

t(t − 1)

=

1 · 2

1 · (−1)

2 · 1

=

y0

(t2 − 3t + 2) − y1(t2

− 2t) +

y2

(t2 − t),

2

2

64

где t = (x − x0)/h. Для оценки погрешности интерполяции получается

выражение

r3(x) = f′′′(ξ(x)) t(t − 1)(t − 2)h3, 3!

где x0 < ξ(x) < x2.

Формула для приближенного вычисления первой производной примет

вид

f′(x) =

f′(x

0

+ th)

1 2t

−

3

2t

−

1

t

y

(2t

2) y

+

y

,

h

≈ h (

2

−

−

2

0

1

2)

а оценка ее погрешности

r3′

(r3(x0 + th))′

f(4)(ξ(x0 + th))

ξt′(x0 + th) t (t − 1) (t − 2) h2+

(x) =

t

=

h

3!

+

f′′′(ξ(x0

+ th))

((t − 1)(t − 2) + t(t − 1) + t(t − 2)) h2

3!

Из оценки погрешности формулы численного дифференцирования следует, что функция f(x) должна иметь четвертую (на порядок больше) ограниченную производную, а порядок по h стал на единицу меньше по

сравнению с интерполированием.

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

вания гладкости функции f(x) не меняются:

f′(x

)

≈

−3y0 + 4y1 − y2

, (t = 0),

|

r′

(x

)

M3

h2

3

0

2h

3

0

| 6

f′(x

)

≈

−y0 + y2

, (t = 1),

|

r′

(x )

| 6

M3

h2

6

1

2h

3

1

f′(x

)

≈

y0 − 4y1 + 3y2

, (t = 2),

|

r′ (x )

|

6

M3

h2

2h

3

2

3

2

где

M

= max

f′′′(x)

3

x06x6x2 |

| .

65

При n = 3 можно еще вычислить вторую производную

f′′(x) =

f′′(x

0

+ th) y

0

−

2y

1

+ y

2

tt

≈

,

h2

h2

оценка погрешности этой формулы имеет вид

((r3(x0 + th))′)′

r3′′(x) =

t

t

= f′′′(ξ(x0 + th))(t − 1)h+

h

+f(5)(ξ(x0 + th))

(t3 − 3t2 + 2)

(ξ′

(x0 + th))2h+

(ξ(x0 + th)) (

3!

t

+ th))h.

t3

3t2 + 2t

(3t2

6t + 2)

+f(4)

−

ξtt′′(x0 + th) +

−

ξt′(x0

3!

3!

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

5.2Метод неопределенных коэффициентов

Метод неопределенных коэффициентов приводит к тем же формулам и имеет два варианта реализации. Этот метод удобен для неравномерно расположенных узлов, вообще в случаях, когда не просто выписывается интерполяционный многочлен. Формулу для вычисления k-ой производной запишем в виде

n
∑i		(5.1)
f(k)(x) ≈	Cif(xi),
=1

здесь xi, (i = 0, 1, ..., n −1) заданные узлы, а Ci пока неизвестные коэффициенты.

Первый вариант метода

Выберем коэффициенты Ci так, чтобы формула (5.1) для вычисления k-ой производной была точна для произвольного многочлена Pm(x) как можно более высокой степени (максимальная алгебраическая точность). Это требование выполняется, если формула точна для любой степенной

66

функции xj (j = 0, 1, ..., m), то есть Ci должны удовлетворять системе линейных алгебраических уравнений (СЛАУ)

n
∑i
(xj)(k) = Cixij, j = 0, 1, ..., m.	(5.2)
=0

Если m = n − 1, то получается неоднородная СЛАУ, состоящая из n уравнений с n неизвестными, и с определителем (Вандермонда ) отличным от нуля. Единственное решение Ci, i = 1, ..., n такой системы определяет формулу приближенного вычисления f(k)(x).

Пример

Пусть n = 3, а узлы x0 = −h, x1 = 0, x2 = h, тогда система (5.2) для

построения формулы вычисления первой производной примет вид

0 = C1

2

+C2

+C3,

2

1 = − C1 h +C2 0 +C3h,

2x = C1 h +C2 0 +C3h .

Вычислив

решения этой системы, получим формулу для приближенного

вычисления первой поизводной

f′(x) =

2x + h

y

2x

y

+

2x − h

y

2 h2

1 − h2

2

2 h2

3

Второй вариант метода.

Каждую функцию f(xi) из выражения (5.1) представим в виде

f(xi) = f(x + (xi − x)) i = 1, 2, ..., n,

затем разложим в ряд Тейлора в окрестности точки x. Подставим эти разложения в (5.1), приведем подобные, приравняем коэффициенты при производных f(m)(x), m ≠ k нулю, а при производной f(k)(x) единице. Получим СЛАУ для вычисления коэффициентов Ci. Число неизвестных уравнений этой системы n не должно превосходить число уравнений m 6 n. Если порядок производной совпадает с числом разыскиваемых коэффициентов

67

f′′(xi) ≈

Ci (k = n), то погрешность построенной формулы будет порядка O(hk+1),

где h = max |x − xi|. Если k < n, то можно получить порядок O(hn+1).

16i6n

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

Пример

Пусть требуется вычислить коэффициенты в выражении

C1 f(xi−1) + C2 f(xi+1) + C3 f(xi+2), (5.3) h2

так, чтобы погрешность этой формулы имела наибольший возможный порядок O(hs) . Здесь узлы считаются равноотстоящими xi − xi−1 = h.

Разложим каждую функцию из этой формулы в окрестности точки xi

в ряд Тейлора

f(xi±1) = f(xi ± h) = f(xi) ±

f′(xi)

f′′(xi)

2

±

f′′′(xi)

3

h +

h

h

+ ... ,

1!

2!

3!

f′(xi)

f′′(xi)

f′′′(xi)

f(xi+2) = f(xi + 2h) = f(xi) +

2h +

4h2

+

8h3 + ... ,

1!

2!

3!

Подставим эти разложения в выражение (5.3) для производной, получим

C1 f(xi−1) + C2 f(xi+1) + C3 f(xi+2)				=	C1 + C2 + C3		f(xi)+
	h2				h2
+	−C1 + C2 + 2C3	f′(xi) +	C1 + C2 + 4C3			f′′(xi)+
	h		2

+−C1 + C2 + 8C3 hf′′′(xi) + . . . .

6

Если выбрать константы так, что

C1 + C2 + C3 = 0

−C1 + C2 + 2C3 = 0

C1 + C2 + 4C3 = 2

то у формулы (5.3) для вычисления второй производной будет погрешность

O(h), так как −C1 + C2 + 8C3 ≠ 0.

68

5.3О вычислительной погрешности численного дифференцирования

Если функции близки, то это не означает, что близки и их производные. Геометрически это видно из рисунка

Рассмотрим влияние вычислительной погрешности (ошибки округления) при численном дифференцировании на примере простой формулы

f′(x0) ≈ y1 − y0 h

Сначала вычислим погрешность самой формулы (погрешность метода). Для этого используем разложение в ряд Тейлора

f(x1) = y1 = y(x0 + h) = y(x0) + hy′(x0) + y′′(ξ)

h2

,

2

где x0 < ξ < x0 + h. Отсюда получается

f′(x

) =

y1 − y0

−

y′′(ξ)

h

,

0

h

2

т.е. рассматриваемая формула имеет погрешность

R

y′′ ξ h

M2h,

M

= max

y′′(x) .

|

1| =

( )

2

6

2

2

x06x6x0+h |

|

Эта оценка погрешности

верна

при отсутствии ошибок округления. Пусть

значения y0 и y1 заданы неточно, а с ошибками δ0 и δ1 причем |δ0| 6 E и

|δ1| 6 E

(y0)T = y0 + δ0, (y1)T = y1 + δ1,

69

причем δ0 6 E и δ1 6 E.

Тогда к ошибке метода добавляется вычислительная погрешность. Вы-

числим ее

(y1)T − (y0)T

=

y1 + δ1 − y0 − δ0

=

y1 − y0

+

δ1 − δ0

.

h

h

h

h

Отсюда получается оценка вычислительной погрешности

| 2|

h

6

h

R

=

δ1 −

δ0

2E .

Теперь общая погрешность при

вычислении

производной

|R| 6 |R1| + |R2| 6

M2h

2E

R = R1 + R2

+

.

2

h

Оценка погрешности Φ(h) =

M2h

2E

имеет минимум Φ(h) > √

,

+

M2 E

2

h

который достигается при h = 2√

E

, то есть выбрать шаг, который га-

M2

рантировал бы погрешность численного дифференцирования меньше чем

√

2 M2 E, нельзя.

При записи чисел в память на ЭВМ относительная ошибка округления имеет оценку E 6 β1−t/2 = 2−t (t– точность ЭВМ, β = 2), тогда

гарантированная оценка погрешности при численном дифференцировании

|R| 6 2−t=2, то есть может потеряться половина верных знаков. Для про-

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

Список литературы

[1]Бахвалов H.С. Численные методы.М.:Hаука,1973.

[2]Бахвалов H.С., Жидков Н.П., Кобельков Г.М. Численные методы. М.:Hаука,1989г.

[3]Форсайт Дж., Малькольм М., Моулер К. Машинные методы математических вычислений. - М.: Мир, 1960.

[4]Калиткин H.H. Численные методы.М.:Hаука,1978

[5]Самарский А.А., Гулин А.В. Численные методы. М.:Hаука,1989.

[6]Вержбицкий В.М. Основы численных методов.- М.:Высш.шк.,2002.

70