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.