Choliy.NumericaMethods.2011Dec01

Метод невизначених коефiцiєнтiв. Формули для чисельного диференцiювання можна отримати також методом невизначених коефiцiєнтiв. Для зручностi вважаємо, що вузли рiвновiддаленi, а точка, в якiй шукається похiдна позначена як x0.

Приймемо для k-ї похiдної такий вираз:

i будемо вимагати, щоб ця формула була точною для полiномiв найвищого можливого степеня. Замiнюємо вираз для функцiї на полiном:

m

X

f(x) = ajxj,

j=0

продиференцiюємо його k разiв:

Якщо тепер прирiвняти доданки з однаковими j, маємо систему (n + 1)-го лiнiйного рiвняння:

n

j(j − 1)(j − 2) . . . (j − k + 1)xj0−k = Xcixji ,

i=1

101

з (n + 1) невiдомими. Розв’язання цiєї системи дасть коефiцiєнти ci для використання у формулi чисельного диференцiювання. Це можна зробити при m + 1 = n, тобто, якщо в таблицi n вузлiв, то це дає змогу побудувати формулу для диференцiювання точну для полiномiв степеня не вище n − 1.

Приклад. Побудуємо формулу для диференцiювання у виглядi: f0(0) = c1f(−h) + c2f(0) + c3f(h).

Для цього випадку система рiвнянь

j(j − 1)(j − 2) . . . (j − k + 1)xj0−k = c1f(−h) + c2f(0) + c3f(h).

У цiй системi k = 1 а у добутку в лiвiй частинi залишається лише один множник, j:

розв’язок якої: c1 = −1/(2h), c2 = 0, c3 = 1/(2h) дозволяє записати досить

вiдому формулу:

f0(0) = f(h) − f(−h). 2h

Приклад. Спробуємо знайти вираз для другої похiдної у такому ж виглядi: f00(0) = c1f(−h) + c2f(0) + c3f(h).

Тепер k = 2 i лiвий добуток тепер мiстить два множники, j(j − 1):

що дає: c1 = c3 = 1/h2, c2 = −2/h2. Формула для диференцiювання така:

f00(0) = f(−h) − 2f(0) + f(h). h2

На цьому прикладi оцiнимо точнiсть формули. Для цього розкладем f(x) в околицi точки x = 0 i пiдставимо цi розвиття у формулу:

R(h) = f00(0) − f(−h) − 2f(0) + f(h) = h2

= f00(0) − f(0) + f0(0)(−h) + f00(0)h2/2 + f000(ξ−)(−h)3/6 − 2f(0)+ +f(0) + f0(0)h + f00(0)h2/2 + f000(ξ+)h3/6 /h2 =

= (f000(ξ−)(−h)3/6 + f000(ξ+)h3/6)/h2 = f(iv)(ξ)h2/3.

Порядок отриманої формули h2.

102

Задачi. У всiх задачах цього роздiлу потрiбно побудувати таблицю похiдних заданої функцiї на iнтервалi вiд 0 до 1 з кроком 0.1 i порiвняти отриманi значення з похiдними, отриманими шляхом диференцiювання заданої формули. Якщо в умовi задачi є параметр, вiн може приймати довiльне значення.

103

12. Iнтегрування

У цьому роздiлi розглядається задача наближеного знаходження визначених iнтегралiв:

b

Z

f(x)dx = F (b) − F (a),

a

звичайно, функцiя f(x) повинна бути на [a, b] неперервною. Чисельно iнтеграли знаходимо коли немає первiсної, або якщо f(x) задана таблицею. Формули, що при цьому використовуються, називаються квадратурними формулами. Як i у попередньому роздiлi, пiдiнтегральну функцiю замiняємо її iнтерполяцiйним полiномом i iнтегруємо його.

Загальна квадратурна формула. Представимо пiдiнтегральну функцiю iнтерполяцiйним полiномом:

ZZ

Подiлимо iнтервал [a, b] на n пiдiнтервалiв: [x0, xn], i = 0, . . . , n−1, i замiнимо пiдiнтегральну функцiю полiномом Лагранжа на такiй сiтцi iнтерполяцiї:

f(x) = Ln(x) + похибка,

тому:

ZZ

f(x)dx = Ln(x)dx + Rn(x),

Rn(x) - iнтегральна похибка. Якщо знехтувати похибкою:

104

при цьому коефiцiєнти

не залежать вiд x, однак, залежать вiд xi, якi виступають як сталi. Таким способом ми отримали загальний вигляд квадратурних формул:

Оскiльки Ai є полiномом, то (12.1) буде точною для полiномiв, степеня, не вищого n. Взявши за пiдiнтегральну функцiю послiдовно полiноми степенiв вiд 0 до n:

побудуємо систему n + 1 лiнiнйних рiвнянь з n + 1 невiдомими Ai, пiсля розв’язання якої квадратурна формула (12.1) стає повнiстю визначеною.

Приклад. Нехай

1

Z

y(x)dx = A0y(1/4) + A1y(1/2) + A2y(3/4).

0

Невiдомих - три коефiцiєнти, тому потрiбно три рiвняння, якi можна побудувати, взявши три перших полiноми:

A0 = 2/3, A1 = −1/3, A2 = 2/3,

105

тобто

Формули Ньютона-Котеса. Це квадратурнi формули для рiвновiддалених вузлiв. У попередньому роздiлi ми отримали вираз для полiному Лагранжа для такого випадку:

Формула для знаходження коефiцiєнтiв така:

0

А загальна квадратурна формула приймає вигляду:

З лiтератури вiдомо досить багато спрощених варiантiв цiєї загальної формули. Отримаємо деякi з них.

Формула трапецiй. Отримаємо, поклавши n = 1:

106

Покажемо на цьому прикладi як знайти похибку iнтегрування. Отже:

Чи трохи по-iншому:

x0+h

Z

h

R(h) = y(x)dx − 2 (y(x0 + h) + y(x0)).

x0

Вiзьмемо похiднi по h:

R0(h) = y(x0+h)−12(y(x0+h)+y(x0))−h2 y0(x0+h) = 12y(x0+h)−12y(x0)−h2 y0(x0+h).

R00(h) = 12y0(x0 + h) − 12y0(x0 + h) − h2 y00(x0 + h) = −h2 y00(x0 + h).

А тепер проiнтегруємо:

де ξ - деяке значення всерединi iнтервалу. Як бачимо, квадратурна формула має порядок h3.

Формулу Сiмпсона можна отримати при n = 2:

107

тобто

x2

Z

h

y(x)dx = 3 (y0 + 4y1 + y2) .

x0

Залишковий член (похибку) знайдемо аналогiчно:

x2

Z

h

R(h) = y(x)dx − 3 (y0 + 4y1 + y2) =

Хоча порядок формули (n = 3) бiльший, нiж формули Сiмпсона (n = 2), остання формула ”гiрша“, оскiльки має бiльшу похибку: 3/80 > 1/90.

Якщо достатньо, щоб квадратурна формула була точна лише для полiномiв нульового степеня, то можна користуватися формулою прямокутникiв:

x1

Z

y(x)dx = hy(x0).

x0

Можна отримати i формули вищих порядкiв, наприклад:

R

y(x)dx = (7y0 + 32y1 + 12y2 + 32y3 + 7y4) /90,

R

y(x)dx = (41y0 + 216y1 + 27y2 + 272y3 + 27y4 + 216y5 + 41y6) /840.

108

Формули Чебишева. Задачу iнтегрування можна поставити дещо по-iншо- му. На вiдмiну вiд формул Ньютона-Котеса, коли вузли сiтки були фiксованi, а квадратурнi формули будувалися шляхом знаходження коефiцiєнтiв, при iншiй постановцi задачi коефiцiєнти вважають вiдомими, а положення вузлiв - шуканими:

Нехай усi Bi однаковi, тодi якщо взяти f(x) = 1,

Формули Гауса. У формулах Ньютона-Котеса використовують змiннi ваги, у формулах Чебишева - пiдбиралися оптимальнi, (змiннi) вузли. При побудовi формул Гауса - i ваги, i вузли вважаються вiльними i їх положення та значення потрiбно визначати, тобто потрiбно так пiдiбрати вузли i ваги, щоб рiвняння

було точним для полiномiв найбiльшого можливого степеня. Формула мiстить Ai, ti - всього 2n невiдомих, значить найбiльша степiнь полiнома, для якого квадратурна формула має бути точною - 2n − 1. Нехай це виконується для

109

1, t, t2, . . . , t2n−1, тодi:

Будь-яку iншу функцiю будемо представляти так:

2n−1

X

f(t) = Cktk,

k=0

Останнє спiввiдношення, переписане у виглядi

є системою нелiнiйних рiвнянь, розв’язання якої дозволяє побудувати квадратурну формулу Гауса. Ця система досить складна, а її розв’язання зовсiм нетривiальне. Спробуємо спростити цей процес. Як i при отриманнi формул Чебишева:

Покладемо тепер

f(x) = tkPn(t), k [0, n − 1],

де Pn(t) - деякий полiном степеня n. Оскiльки f(x) є полiномом, степеня, не вище 2n − 1, то f(x) повиннi задовольняти спiввiдношення (12.2):

110