8: Інтерполяція

Розглянемо наступну задачу: За даними значень відомої функції\(y=f(x)\) в послідовності впорядкованих точок\(x_{0}, x_{1}, \ldots, x_{n}\) знайдіть\(f(x)\) для довільного\(x .\) Коли\(x_{0} \leq x \leq x_{n}\), задача називається інтерполяцією. Коли\(x<x_{0}\) або\(x>x_{n}\), проблема називається екстраполяцією.

З\(y_{i}=f\left(x_{i}\right)\), проблема інтерполяції в основному полягає в тому, щоб промальовувати плавну криву через відомі точки\(\left(x_{0}, y_{0}\right),\left(x_{1}, y_{1}\right), \ldots,\left(x_{n}, y_{n}\right) .\) Це не та сама проблема, що малювання гладкої кривої, яка наближає набір точок даних, які мають експериментальну помилку. Ця остання задача називається наближенням найменших квадратів, яка розглядається в наступному розділі.

Можна інтерполювати\(n+1\) відомі точки унікальним поліномом ступеня\(n\). Коли\(n=1\), многочлен є лінійною функцією; коли\(n=2\), многочлен є квадратичною функцією. Хоча поліноми низького порядку іноді використовуються, коли кількість точок мало, інтерполяції поліномів вищого порядку, як правило, нестабільні і не мають особливого практичного використання.

Тут ми розглянемо більш корисну кусково-поліноміальну інтерполяцію. Двома найбільш використовуваними є кусково-лінійна інтерполяція та кубічна інтерполяція сплайнів. Перший використовує лінійні многочлени, а другий кубічні поліноми.

Кусково-лінійна інтерполяція

Тут ми використовуємо лінійні многочлени. Це типова інтерполяція, яка зазвичай використовується при побудові даних.

Припустимо\(y=g(x)\), що функція інтерполяції є, і, як і раніше, є\(n+1\) точки для інтерполяції. Побудовано функцію\(g(x)\) з\(n\) локальних лінійних многочленів. пишемо

\[g(x)=g_{i}(x), \quad \text { for } x_{i} \leq x \leq x_{i+1}, \nonumber \]

де

\[g_{i}(x)=a_{i}\left(x-x_{i}\right)+b_{i} \nonumber \]

і\(i=0,1, \ldots, n-1\)

Тепер ми\(y=g_{i}(x)\) вимагаємо пройти через кінцеві точки\(\left(x_{i}, y_{i}\right)\) і\(\left(x_{i+1}, y_{i+1}\right) .\) Ми маємо

\[\begin{aligned} y_{i} &=b_{i} \\ y_{i+1} &=a_{i}\left(x_{i+1}-x_{i}\right)+b_{i} \end{aligned} \nonumber \]

Тому коефіцієнти\(g_{i}(x)\) визначаються як

\[a_{i}=\frac{y_{i+1}-y_{i}}{x_{i+1}-x_{i}}, \quad b_{i}=y_{i} \nonumber \]

Хоча кусково-лінійна інтерполяція широко використовується, особливо в побудові процедур, вона страждає від розриву похідної в кожній точці. Це призводить до функції, яка може виглядати не гладкою, якщо точки занадто широко розташовані. Далі ми розглянемо більш складний алгоритм, який використовує кубічні многочлени.

Кубічна інтерполяція сплайнів

\(n+1\)Точки, які потрібно інтерполювати, знову

\[\left(x_{0}, y_{0}\right),\left(x_{1}, y_{1}\right), \ldots\left(x_{n}, y_{n}\right) . \nonumber \]

Тут ми використовуємо\(n\) кусково-кубічні многочлени для інтерполяції,

\[g_{i}(x)=a_{i}\left(x-x_{i}\right)^{3}+b_{i}\left(x-x_{i}\right)^{2}+c_{i}\left(x-x_{i}\right)+d_{i}, \quad i=0,1, \ldots, n-1, \nonumber \]

з функцією глобальної інтерполяції, записаною як

\[g(x)=g_{i}(x), \quad \text { for } x_{i} \leq x \leq x_{i+1} . \nonumber \]

Для досягнення плавної інтерполяції ми накладаємо, що\(g(x)\) і його перша та друга похідні є безперервними. Вимога, яка\(g(x)\) є безперервною (і проходить через всі\(n+1\) пункти) призводить до двох обмежень

\[\begin{align} g_{i}\left(x_{i}\right) &=y_{i}, \quad i=0 \text { to } n-1 \\ g_{i}\left(x_{i+1}\right) &=y_{i+1}, \quad i=0 \text { to } n-1 \end{align} \nonumber \]

Вимога, яка\(g^{\prime}(x)\) є безперервним, призводить до

\[g_{i}^{\prime}\left(x_{i+1}\right)=g_{i+1}^{\prime}\left(x_{i+1}\right), \quad i=0 \text { to } n-2 \nonumber \]

І вимога, яка\(g^{\prime \prime}(x)\) є постійною, призводить до

\[g_{i}^{\prime \prime}\left(x_{i+1}\right)=g_{i+1}^{\prime \prime}\left(x_{i+1}\right), \quad i=0 \text { to } n-2 . \nonumber \]

Існують\(n\) кубічні многочлени,\(g_{i}(x)\) і кожен кубічний многочлен має чотири вільних коефіцієнта; тому існує загальна кількість\(4 n\) невідомих коефіцієнтів. Число стримуючих рівнянь з (8.7) - (8.10) дорівнює\(2 n+2(n-1)=4 n-2\). При\(4 n-2\)\(4 n\) обмеженнях і невідомих для унікального рішення потрібні ще дві умови. Зазвичай вони вибираються як додаткові умови на першому\(g_{0}(x)\) та останньому\(g_{n-1}(x)\) многочленах. Ці додаткові умови ми обговоримо пізніше.

Тепер перейдемо до визначення рівнянь для невідомих коефіцієнтів кубічних многочленів. Поліноми та їх перші дві похідні задаються

\[\begin{align} g_{i}(x) &=a_{i}\left(x-x_{i}\right)^{3}+b_{i}\left(x-x_{i}\right)^{2}+c_{i}\left(x-x_{i}\right)+d_{i} \\ g_{i}^{\prime}(x) &=3 a_{i}\left(x-x_{i}\right)^{2}+2 b_{i}\left(x-x_{i}\right)+c_{i} \\ g_{i}^{\prime \prime}(x) &=6 a_{i}\left(x-x_{i}\right)+2 b_{i} \end{align} \nonumber \]

Ми розглянемо чотири умови (8.7) - (8.10) по черзі. З (8.7) і (8.11) ми маємо

\[d_{i}=y_{i}, \quad i=0 \text { to } n-1, \nonumber \]

який безпосередньо вирішує для всіх\(d\) -коефіцієнтів.

Щоб задовольнити (8.8), ми спочатку визначаємо

\[h_{i}=x_{i+1}-x_{i}, \nonumber \]

і

\[f_{i}=y_{i+1}-y_{i} . \nonumber \]

Тепер з (8.8) і (8.11), використовуючи (8.14), отримаємо\(n\) рівняння

\[a_{i} h_{i}^{3}+b_{i} h_{i}^{2}+c_{i} h_{i}=f_{i}, \quad i=0 \text { to } n-1 . \nonumber \]

З (8.9) і (8.12) отримуємо\(n-1\) рівняння

\[3 a_{i} h_{i}^{2}+2 b_{i} h_{i}+c_{i}=c_{i+1}, \quad i=0 \text { to } n-2 \nonumber \]

З (8.10) і (8.13) отримуємо\(n-1\) рівняння

\[3 a_{i} h_{i}+b_{i}=b_{i+1} \quad i=0 \text { to } n-2 \text {. } \nonumber \]

Корисно буде включити визначення коефіцієнта\(b_{n}\), якого зараз немає. (Індекс коефіцієнтів кубічних поліномів тільки піднімається до\(n-1\).) Ми просто розширюємо (8.19) до\(i=n-1\) і так пишемо

\[3 a_{n-1} h_{n-1}+b_{n-1}=b_{n}, \nonumber \]

які можна розглядати як визначення\(b_{n}\).

Тепер ми приступимо до усунення множин\(a\) - і c-коефіцієнтів (з\(d\) -коефіцієнтами, які вже усунені в (8.14)), щоб знайти систему лінійних рівнянь для\(b\) -коефіцієнтів. З (8.19) і\((8.20)\), ми можемо вирішити для\(n a\) -коефіцієнтів, щоб знайти

\[a_{i}=\frac{1}{3 h_{i}}\left(b_{i+1}-b_{i}\right), \quad i=0 \text { to } n-1 . \nonumber \]

З (8.17) ми можемо вирішити для\(n\) c-коефіцієнтів наступним чином:

\[\begin{align} \nonumber c_{i} &=\frac{1}{h_{i}}\left(f_{i}-a_{i} h_{i}^{3}-b_{i} h_{i}^{2}\right) \\ \nonumber &=\frac{1}{h_{i}}\left(f_{i}-\frac{1}{3 h_{i}}\left(b_{i+1}-b_{i}\right) h_{i}^{3}-b_{i} h_{i}^{2}\right) \\ &=\frac{f_{i}}{h_{i}}-\frac{1}{3} h_{i}\left(b_{i+1}+2 b_{i}\right), \quad i=0 \text { to } n-1 \end{align} \nonumber \]

Тепер ми можемо знайти рівняння для\(b\) -коефіцієнтів шляхом підстановки\((8.21)\) і\((8.22)\) в\((8.18)\):

\[\begin{array}{r} 3\left(\frac{1}{3 h_{i}}\left(b_{i+1}-b_{i}\right)\right) h_{i}^{2}+2 b_{i} h_{i}+\left(\frac{f_{i}}{h_{i}}-\frac{1}{3} h_{i}\left(b_{i+1}+2 b_{i}\right)\right) \\ =\left(\frac{f_{i+1}}{h_{i+1}}-\frac{1}{3} h_{i+1}\left(b_{i+2}+2 b_{i+1}\right)\right), \nonumber \end{array} \nonumber \]

що спрощує

\[\frac{1}{3} h_{i} b_{i}+\frac{2}{3}\left(h_{i}+h_{i+1}\right) b_{i+1}+\frac{1}{3} h_{i+1} b_{i+2}=\frac{f_{i+1}}{h_{i+1}}-\frac{f_{i}}{h_{i}} \nonumber \]

рівняння, яке є дійсним\(i=0\) для\(n-2\). Тому (8.23) представляють\(n-1\) рівняння для\(n+1\) невідомих\(b\) -коефіцієнтів. Відповідно, запишемо матричне рівняння для\(b\) -коефіцієнтів, залишивши перший і останній рядок відсутніми, як

\[\left(\begin{array}{cccccccc} \ldots & \ldots & \ldots & \ldots & \text { missing } & \ldots & \ldots \\ \frac{1}{3} h_{0} & \frac{2}{3}\left(h_{0}+h_{1}\right) & \frac{1}{3} h_{1} & \ldots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \ldots & \frac{1}{3} h_{n-2} & \frac{2}{3}\left(h_{n-2}+h_{n-1}\right) & \frac{1}{3} h_{n-1} \\ \ldots & \ldots & \ldots & \ldots & \text { missing } & \ldots & \ldots \end{array}\right)\left(\begin{array}{c} b_{0} \\ b_{1} \\ \vdots \\ b_{n-1} \\ b_{n} \end{array}\right)\left(\begin{array}{c} \text { missing } \\ \frac{f_{1}}{h_{1}}-\frac{f_{0}}{h_{0}} \\ \vdots \\ \frac{f_{n-1}}{h_{n-1}}-\frac{f_{n-2}}{h_{n-2}} \\ \text { missing } \end{array}\right) . \nonumber \]

Після того, як задані відсутні перше та останнє рівняння, матричне рівняння для\(b\) -коефіцієнтів може бути вирішено шляхом гаусової елімінації. І як тільки\(b\) -коефіцієнти визначені,\(a\) - і\(c\) -коефіцієнти також можуть бути визначені з (8.21) і (8.22), з\(d\) -коефіцієнтами, вже відомими з (8.14). Кусково-кубічні многочлени, отже, відомі і\(g(x)\) можуть бути використані для інтерполяції до будь-якого\(x\) задовольняє значення\(x_{0} \leq x \leq x_{n}\).

Відсутні перше і останнє рівняння можуть бути вказані декількома способами, і тут ми покажемо два способи, які дозволені функцією MATLAB spline.m. Перший спосіб слід використовувати, коли похідна\(g^{\prime}(x)\) відома в кінцевих точках,\(x_{0}\) і\(x_{n} ;\) тобто припустимо, ми знаємо значення\(\alpha\) і\(\beta\) такі, що

\[g_{0}^{\prime}\left(x_{0}\right)=\alpha, \quad g_{n-1}^{\prime}\left(x_{n}\right)=\beta . \nonumber \]

З відомого значення\(\alpha\) і використання (8.12) і (8.22), ми маємо

\[\begin{aligned} \alpha &=c_{0} \\ &=\frac{f_{0}}{h_{0}}-\frac{1}{3} h_{0}\left(b_{1}+2 b_{0}\right) \end{aligned} \nonumber \]

Тому відсутнє перше рівняння визначається як

\[\frac{2}{3} h_{0} b_{0}+\frac{1}{3} h_{0} b_{1}=\frac{f_{0}}{h_{0}}-\alpha \nonumber \]

Від відомого значення\(\beta\), і використання\((8.12),(8.21)\), і\((8.22)\), ми маємо

\[\begin{aligned} \beta &=3 a_{n-1} h_{n-1}^{2}+2 b_{n-1} h_{n-1}+c_{n-1} \\ &=3\left(\frac{1}{3 h_{n-1}}\left(b_{n}-b_{n-1}\right)\right) h_{n-1}^{2}+2 b_{n-1} h_{n-1}+\left(\frac{f_{n-1}}{h_{n-1}}-\frac{1}{3} h_{n-1}\left(b_{n}+2 b_{n-1}\right)\right), \end{aligned} \nonumber \]

що спрощує

\[\frac{1}{3} h_{n-1} b_{n-1}+\frac{2}{3} h_{n-1} b_{n}=\beta-\frac{f_{n-1}}{h_{n-1}} \nonumber \]

буде використано як відсутнє останнє рівняння.

Другий спосіб вказівки відсутніх першого і останнього рівнянь називається умовою not-a-knot, яке передбачає, що

\[g_{0}(x)=g_{1}(x), \quad g_{n-2}(x)=g_{n-1}(x) \nonumber \]

Розглядаючи перше з цих рівнянь, з (8.11) ми маємо

\[a_{0}\left(x-x_{0}\right)^{3}+b_{0}\left(x-x_{0}\right)^{2}+c_{0}\left(x-x_{0}\right)+d_{0} \nonumber \]

\[=a_{1}\left(x-x_{1}\right)^{3}+b_{1}\left(x-x_{1}\right)^{2}+c_{1}\left(x-x_{1}\right)+d_{1} . \nonumber \]

Тепер два кубічні многочлени можуть бути ідентичні\(x\), якщо при певному значенні поліноми та три їх перші похідні ідентичні. Наші умови безперервності\(x=x_{1}\) вже вимагають, щоб при цьому значенні\(x\) цих двох поліномів і їх перших двох похідних були однаковими. Самі многочлени будуть ідентичними, то, якщо їх треті похідні теж ідентичні при\(x=x_{1}\), або якщо

\[a_{0}=a_{1} . \nonumber \]

З (8.21) ми маємо

\[\frac{1}{3 h_{0}}\left(b_{1}-b_{0}\right)=\frac{1}{3 h_{1}}\left(b_{2}-b_{1}\right), \nonumber \]

або після спрощення

\[h_{1} b_{0}-\left(h_{0}+h_{1}\right) b_{1}+h_{0} b_{2}=0, \nonumber \]

який надає нам наше відсутнє перше рівняння. Подібний аргумент\(x=x_{n-1}\) також надає нам наше останнє рівняння,

\[h_{n-1} b_{n-2}-\left(h_{n-2}+h_{n-1}\right) b_{n-1}+h_{n-2} b_{n}=0 . \nonumber \]

Підпрограми MATLAB spline.m і ppval.m можуть бути використані для кубічної інтерполяції сплайнів (див. також interp1.m). Я проілюструю ці процедури в класі і розмістити зразок коду на веб-сайті курсу.

Багатовимірна інтерполяція

Припустимо, ми інтерполяція значення функції двох змінних,

\[z=f(x, y) . \nonumber \]

Відомі значення задаються

\[z_{i j}=f\left(x_{i}, y_{j}\right), \nonumber \]

з\(i=0,1, \ldots, n\) і\(j=0,1, \ldots, n .\) Зверніть увагу, що\((x, y)\) точки, в яких\(f(x, y)\) відомі, лежать на сітці в\(x-y\) площині.

\(z=g(x, y)\)Дозволяти інтерполяція функція,\(z_{i j}=g\left(x_{i}, y_{j}\right) .\) задовольняє двовимірну інтерполяцію, щоб знайти значення\(g\) в точці\((x, y)\) може бути зроблено шляхом першого виконання\(n+1\) одновимірних інтерполяцій в \(y\)знайти значення\(g\) в\(n+1\) точках\(\left(x_{0}, y\right),\left(x_{1}, y\right), \ldots\)\(\left(x_{n}, y\right)\), з подальшою одновимірною інтерполяцією\(x\) в знайти значення\(g\) at\((x, y)\).

Іншими словами, двовимірна інтерполяція на розмірній сітці\((n+1) \times(n+1)\) виконується спочатку\(n+1\) одновимірними інтерполяціями до значення з\(y\) подальшою одновимірною інтерполяцією до значення\(x\). Двовимірна інтерполяція може бути узагальнена до більш високих вимірів. Функціями MATLAB для виконання дво- і тривимірної інтерполяції є interp2.m і interp3.m.