Ошибка интерполяции кубическими сплайнами

Материал из MachineLearning.

Перейти к: навигация, поиск

Содержание

1 Введение
- 1.1 Постановка математической задачи
2 Изложение метода
- 2.1 Метод прогонки
- 2.2 Пример: интерполирование неизвестной функции
3 Ошибка интерполяции
- 3.1 Пример: интерполяция синуса
4 Список литературы
5 См. также

Введение

Постановка математической задачи

Одной из основных задач численного анализа является задача об интерполяции функций.
Пусть на отрезке $ale xle b$ задана сетка $omega={x_i:x_0=a<x_1<cdots<x_i<cdots<x_n=b}$ и в её узлах заданы значения функции $y(x)$ , равные $y(x_0)=y_0,cdots,y(x_i)=y_i,cdots,y(x_n)=y_n$ . Требуется построить интерполянту — функцию $f(x)$ , совпадающую с функцией $y(x)$ в узлах сетки:

( 1 )

$f(x_i)=y_i, i=0,1,cdots,n.$

Основная цель интерполяции — получить быстрый (экономичный) алгоритм вычисления значений $f(x)$ для значений $x$ , не содержащихся в таблице данных.

Интерполируюшие функции $f(x)$ , как правило строятся в виде линейных комбинаций некоторых элементарных функций:

$f(x)=sum_{k=0}^N {c_kPhi_k(x)},$

где ${Phi_k(x)}$ — фиксированный линейно независимые функции, $c_0, c_1, cdots, c_n$ — не определенные пока коэффициенты.

Из условия (1) получаем систему из $n+1$ уравнений относительно коэффициентов ${c_k}$ :

$sum_{k=0}^N {c_kPhi_k(x_i)}=y_i, i=0,1,cdots,n.$

Предположим, что система функций $Phi_k(x)$ такова, что при любом выборе узлов $a<x_1<cdots<x_i<cdots<x_n=b$ отличен от нуля определитель системы:

$Delta(Phi) = begin{vmatrix} Phi_0(x_0) & Phi_1(x_0) & cdots & Phi_n(x_0) \Phi_0(x_1) & Phi_1(x_1) & cdots & Phi_n(x_1)\ cdots & cdots & cdots & cdots \Phi_0(x_n) & Phi_1(x_n) & cdots & Phi_n(x_n)end{vmatrix}$ .

Тогда по заданным $y_i (i=1,cdots, n)$ однозначно определяются коэффициенты $c_k (k=1,cdots, n)$ .

Изложение метода

Интерполяция кубическими сплайнами является частным случаем кусочно-полиномиальной интерполцией. В этом специальном случае между любыми двумя соседними узлами функция интерполируется кубическим полиномом. его коэффициенты на каждом интервале определяются из условий сопряжения в узлах:

$f_i=y_i, f'(x_i-0)=f'(x_i+0), f''(x_i-0)=f''(x_i+0), i=1, 2, cdots, n-1.$

Кроме того, на границе при $x=x_0$ и $x=x_n$ ставятся условия

( 2 )

$f''(x_0)=0, f''(x_n)=0.$

Будем искать кубический полином в виде

( 3 )

$f(x)=a_i+b_i(x-x_{i-1})+c_i(x-x_{i-1})^2+d_i(x-x_{i-1})^3, x_{i-1}le xle x_i.$

Из условия $f_i=y_i$ имеем

( 4 )

$f(x_{i-1})=a_i=y_{i-1}, f(x_i)=a_i+b_ih_i+c_ih_i^2+d_ih_i^3=y_i, h_i=x_i-x_{i-1}, i=1, 2, cdots, n-1.$

Вычислим производные:

$f'(x)=b_i+2c_i(x-x_{i-1})+3d_i(x-x_{i-1})^2, f''(x)=2c_i+6d_i(x-x_{i-1}), x_{i-1}le xle x_i,$

и потребуем их непрерывности при $x=x_i$ :

( 5 )

$b_{i+1}=b_i+2c_ih_i + 3d_ih_i^2, c_{i+1}=c_i+3d_ih_i, i=1, 2, cdots, n-1.$

Общее число неизвестных коэффициентов, очевидно, равно $4n$ , число уравнений (4) и (5) равно $4n-2$ . Недостающие два уравнения получаем из условия (2) при $x=x_0$ и $x=x_n$ :

$c_1=0, c_n+3d_nh_n=0.$

Выражение из (5) $d_i=frac{c_{i+1}-c_i}{3h_i}$ , подставляя это выражение в (4) и исключая $a_i=y_{i-1}$ , получим

$b_i=[frac{y_i-y_{i-1}}h_i]-frac{1}{3}h_i(c_{i+1}+2c_i), i=1, 2, cdots, n-1, b_n=[frac{y_n-y_{n-1}}h_n]-frac{2}{3}h_nc_n,.$

Подставив теперь выражения для $b_i, b_{i+1}$ и $d_i$ в первую формулу (5), после несложных преобразований получаем для определения $c_i$ разностное уравнение второго порядка

( 6 )

$h_ic_i+2(h_i+h_{i+1})c_{i+1}+h_{i+1}c_{i+2}=3left(frac{y_{i+1}-y_i}{h_{i+1}} - frac{y_i-y_{i-1}}{h_i}right), i=1, 2, cdots, n-1.$

С краевыми условиями

( 7 )

$c_1=0, c_{n+1}=0.$

Условие $c_{n+1}=0$ эквивалентно условию $c_n+3d_nh_n=0$ и уравнению $c_{i+1} = c_i+d_ih_i$ . Разностное уравнение (6) с условиями (7) можно решить методом прогонки, представив в виде системы линейных алгебраических уравнений вида $~A*x=F$ , где вектор $x$ соответствует вектору ${c_i}$ , вектор $F$ поэлементно равен правой части уравнения (6), а матрица $~A$ имеет следующий вид:

$A = begin{pmatrix} C_1 & B_1 & 0 & 0 & cdots & 0 & 0 \ A_2 & C_2 & B_2 & 0 & cdots & 0 & 0 \ 0 & A_3 & C_3 & B_3 & cdots & 0 & 0 \ cdots & cdots & cdots & cdots & cdots & cdots & cdots \ cdots & cdots & cdots & cdots & cdots & cdots & cdots \ cdots & cdots & cdots & cdots & cdots & cdots & B_{n-1} \ 0 & 0 & 0 & 0 & cdots & A_{n} & C_{n} end{pmatrix},$

где $A_i=h_i, i=2, cdots, n, B_i = h_{i+1}, i=1, cdots, n-1$ и $C_i=2(h_i+h_{i+1}), i =1, cdots, n$ .

Метод прогонки

Метод прогонки, основан на предположении, что искомые неизвестные связаны рекуррентным соотношением:

( 8 )

$x_i = alpha_{i+1}x_{i+1} + beta_{i+1}, i=1,cdots,n-1$

Используя это соотношение, выразим $x_{i-1}$ и $x_i$ через $x_{i+1}$ и подставим в i-e уравнение:

$left(A_ialpha_ialpha_{i+1} + C_ialpha_{i+1} + B_iright)x_{i+1} + A_ialpha_ibeta_{i+1} + A_ibeta_i + C_ibeta_{i+1} - F_i = 0$

где $F_i$ — правая часть i-го уравнения. Это соотношение будет выполняться независимо от решения, если потребовать

$A_ialpha_ialpha_{i+1} + C_ialpha_{i+1} + B_i = 0$

$A_ialpha_ibeta_{i+1} + A_ibeta_i + C_ibeta_{i+1} - F_i = 0$

Отсюда следует:

$alpha_{i+1} = {-B_i over A_ialpha_i + C_i}$

$beta_{i+1} = {F_i - A_ibeta_i over A_ialpha_i + C_i}$

Из первого уравнения получим:

$alpha_2 = {-B_1 over C_1}$

$beta_2 = {F_1 over C_1}$

После нахождения прогоночных коэффициентов $alpha$ и $beta$ , используя уравнение (1), получим решение системы. При этом,

$x_n = {F_n-A_nbeta_n over C_n+A_nalpha_n}$

Пример: интерполирование неизвестной функции

Построим интерполянту для для функции $f$ , заданной следующим образом:

Вводные значения для задачи интерполяции

$x$	$f(x)$
1	1.0002
2	1.0341
3	0.6
4	0.40105
5	0.1
6	0.23975

В результате интерполяции были рассчитаны следующие коэффициенты интерполянты:

Результат интерполяции

$a$	$b$	$c$	$d$	Отрезок
1,0002	-0,140113846	0,440979231	-0,266965385	$1le xle 2$
1,0341	-0,291901538	-0,359916923	0,217718462	$2le xle 3$
0,6	-0,22553	0,293238462	-0,266658462	$3le xle 4$
0,40105	-0,100328462	-0,506736923	0,306015385	$4le xle 5$
0,1	-0,134456154	0,411309231	-0,137103077	$5le xle 6$

Ошибка интерполяции

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

$parallel s-fparallel = max_{xin[a,b]}|s(x)-f(x)|$

от шага h, где $h = max_{k=1,2,cdots,n-1}|x_{k+1}-x_k|$ .

Известно, что если функция $ƒ(x)$ имеет четыре непрерывные производные, то для ошибки интерполяции определенным выше кубическим сплайном $s(x)$ верна следующая оценка

$parallel s-fparallel le frac{5}{384}h^4parallel frac{d^4f}{df^4}parallel$

причем константа $frac{5}{384}$ в этом неравенстве является наилучшей из возможных

Пример: интерполяция синуса

Постром интерполянту функции $f=sin(4x)$ на отрезке $[-1;1]$ , взяв равномерно отстоящие узлы с шагом 0.5 и шагом 0.25, и сравним полученные результаты.

	Ошибка интерполяции	Оценка ошибки	Иллюстрация
$h=0.5$	0.429685	3.(3)	Результат интерполяции sin(4x) с шагом 0.5
$h=0.25$	0.005167	0.208(3)	Результат интерполяции sin(4x) с шагом 0.25

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

Код программы на языке С++ с помощью которой были произведены все расчеты, можно скачать тут.

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

А.А.Самарский, А.В.Гулин. Численные методы М.: Наука, 1989.
А.А.Самарский. Введение в численные методы М.: Наука, 1982.
Костомаров Д.П., Фаворский А.П. Вводные лекции по численным методам
Н.Н.Калиткин. Численные методы М.: Наука, 1978.

См. также

Интерполяция каноническим полиномом
Тригонометрическая интерполяция
Рациональная интерполяция
Проблема выбора узлов для интерполяции
Практикум ММП ВМК, 4й курс, осень 2008

Источник

bn = 2hn−3 (yn−1 − yn )− 2hn−−31 (yn−2 − yn−1 ). Такие системы быстро решаются специальными методами, с которыми мы познакомимся позднее.

Теорема 3.10. Пусть функция

y = f (x) имеет на отрезке [a, b] непрерывную про-

изводную четвертого порядка и M

= max

f (4)(x)

. Тогда для интерполяционного куби-

[a,b]

ческого сплайна S3 (x), удовлетворяющего граничным условиям типов 1, 2, 4 и 5, спра-

ведлива следующая оценка погрешности:

max

(x)− S

(x)

≤ C M

, где C — некото-

рая константа.

[a, b]

max

не только сам аппроксимирует функцию y = f (x), но

Характерно, что сплайн S3 (x)

и его производные S3/ (x), S3// (x), S3/// (x)

приближают соответствующие производные функции

y = f (x).

Теорема 3.11. При выполнении условий теоремы 3.10 для указанных в ней

сплайнов справедливы неравенства max

(k )(x)− S (k )(x)

≤ C

h4−k , k = 1,2,3.

[a, b]

Пример. Пусть функция задана следующей таблицей:

i		0		1		2		3		4

xi	0		1		2		3		4

yi		1.0		1.8		2.2		1.4		1.0

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

n = 4,

h = 1 = const . Запишем систему (3.10.4), реализующую эту схему в конкретном число-

вом

выражении.

Первое

уравнение

имеет

вид

h1−2 s0 + (h1−2 − h2−2 )s1

− h2−2 s2 =

= 2h−3

− y

)− 2h−3 (y

− y

Вычисляя значения

h − h

−1

и подставляя

h = 1 ,

получим

s0 + 0

s1 − s2 = 2 (− 0.4)− 2 (− 0.8) s0 − s2

= 0.8. Аналогично,

последнее

уравнение дает

h−2 s

n−2

+ (h−2

− h−2 )s

n−1

− h−2 s

= 2h−3

n−1

− y

)− 2h−3

n−2

− y

n−1

), n = 4, n −1 = 3, n − 2 = 3.

n−1

Тогда

h3−2 s2

+ (h3−2

− h4−2 )s3

− h4−2 s4

= 2h4−3 (y3 − y4

)− 2h3−3 (y2

− y3 ) s2 − 0 s3 − s4 =

= 2 0.4 − 2 0.8 s2 − s4 = −0.8. Три внутренние узла дадут три следующих уравнения:

i = 1 : h1−1 s0 + 2 (h1−1 + h2−1 )s1 + h2−1 s 2 = 3[h1−2 (y1 − y 0 )+ h2−2 (y 2 − y1 )],
i = 2 : h2−1 s1 + 2	(h2−1	+ h3−1 )s2	+ h3−1 s3 = 3	[h2−2 (y 2	− y1 )+ h3−2 (y 3 − y 2	)],
i = 3 : h3−1 s 2 + 2	(h3−1	+ h4−1 )s3	+ h4−1 s 4 = 3	[h3−2 (y 3	− y 2 )+ h4−2 (y 4 − y 3	)].

Упрощая их аналогичным способом, что и выше, получим

s0 + 2 2 s1 + s2 = 3(0.8 + 0.4), s1 + 2 2 s2 + s3 = 3(0.4 + (− 0.8)), s2 + 2 2 s3 + s4 = 3(− 0.8 + (− 0.4)).

Тогда вся система имеет вид

− s2

= 0.8,

+ 4s1 + s2

= 3.6,

s1 + 4s2

+ s3

= −1.2,

+ 4s3

+ s4

= −3.6,

−

= −0.8.

Решим эту систему методом Гаусса:

− s2

= 0.8,

− s2

= 0.8,

4s1 + 2s2

= 2.8,

s1 + 4s2

+ s3

= −1.2,

s1 + 4s2

+ s3

= −1.2,

−14s2

− 4s3

= 7.6,

s2 + 4s3 + s

= −3.6,

s2 + 4s3

+ s4 = −3.6,

−

= −0.8.

− s

= −0.8.

− s2

= 0.8,

− s2 = 0.8,

s1 + 4s2

+ s3

= −1.2,

s1 + 4s2 + s3

= −1.2,

+ 4s3 + s4

= −3.6,

s2 + 4s3

+ s4

= −3.6,

− 4s3 − 2s4

= 2.8,

s3 + 0.5s4 = −0.7,

52s

+14s

= −42.8.

−12s

= −6.4.

Отсюда

s4 = 8 /15, s3

= −29 / 30,

= −4 /15, s1

= 5 / 6,

= 8 /15.

Запись формулы (3.10.1)

кубического интерполяционного многочлена Эрмита,

задавая yi−1 ,

yi , si−1 ,

однозначно

определяет сплайн на отрезке [xi−1 , xi ].

Эта формула уже учитывает, что интерполированная

функция

f (x) будет непрерывна и один раз дифференцируема везде на сетке

x0 , x1 ,…, xn .

Таким образом, на каждом из четырех отрезков

[0, 1], [1, 2], [2, 3], [3, 4] конкретный вид

сплайна

можно

установить

по

формуле

(3.10.1).

Например,

для

[0,1]

имеем

y0 = 1.0,

= 1.8,

s0 = 8 / 15,

s1 = 5 / 6. Тогда

(x)

= y

(x1 − x)2 (2(x − x0 )+ h)

+ y /

(x1 − x)2 (x − x0 )

3,1

+ y

(x − x0 )2 (2(x1 − x)

+ h)

+ y /

(x − x0 )2 (x − x1 )

, h = x − x

(x)=1 (1 − x)2 (2(x − 0)+1)

(1 − x)2 (x − 0)

+1.8 (x − 0)2 (2(1 − x)+1)+

3,1

(x − 0)2 (x −1) = (x −1)2 (2x +1)+

(x −1)2 x +1.8x2 (3 − 2x)+

x2 (x −1)=

= −

x3 +

x2 +

x +1 = −0.233333x3 + 0.500000 x2

+ 0.533333x +1.000000.

Многочлены P3,i

на остальных отрезках рассчитываются аналогично. Именно:

x [1, 2],

= 1,

= 2,

= 1.8,

y2 = 2.2,

= 0.833333, s2 = −0.266667,

(x)= −0.233333x3 + 0.500000 x2

+ 0.533333x +1.000000.

3,2

x [2, 3],

= 2, x3

= 3,

= 2.2,

y3 =1.4, s2

= −0.266667, s3

= −0.966667,

(x)= 0.366667 x3 − 3.100000 x2

+ 7.700000 x − 3.800000.

3,3

x [3, 4],

= 3,

= 4,

=1.4,

y4 =1.0,

= −0.966667, s4

= 0.533333,

(x)= 0.366667 x3 − 3.100000 x2

+ 7.700000 x − 3.800000.

3,4

3.12. Лабораторная работа № 6. Аппроксимация функций кубическими сплайнами

Часто на практике, для того чтобы аппроксимировать функцию, вместо построения глобального интерполяционного многочлена на всем промежутке определения аргумента используют интерполяцию кусочными многочленами – сплайнами.

Гладкие сплайны обладают многими замечательными свойствами, которые обеспечили им успех в приложениях. Эти свойства описаны в подразд 3.10. Решим в пакете Mathcad задачу аппроксимации таблично заданной функции с помощью кубических сплайнов.

Среди одиннадцати функций интерполяции и экстраполяции в среде Mathcad имеются средства линейной интерполяции (функция linterp(vx,vy,x) – значение в точке х линейного интерполяционного многочлена векторов vx и vy) и интерполяции сплайнами, осуществляемой функцией interp(K,X,Y,x) по значениям коэффициентов K векторов линейного, параболического или кубического сплайнов. Вектор коэффициентов соответствующего сплайна вычисляется подпрограммами lspline(vx,vy), pspline(vx,vy) и cspline(vx,vy).

Различия в линейной, параболической и кубической интерполяции сплайнами заметно проявляются лишь на концах заданной сетки узлов.

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

	0.593		0.53305
	0.598		0.53464

	0.605		0.54160
ORIGIN := 1	0.613		0.54324
X :=	Y :=	.
	0.619		0.54043
			0.55598
	0.627
	0.632		0.55843

Проведем сначала линейную интерполяцию и построим график результата. y(x):= linterp(X ,Y , x)

Для того чтобы исходные данные представлялись в виде квадратиков, надо войти в панель форматиро-

вания графиков (Formatting Currently Selected X-Y Plot), щелкнув дважды левой кнопкой мыши по любому месту графика. В открывшейся панели

для графика № 2 (trace 2) следует выбрать представление кривой в виде точек (points) квад-

ратной формы (box).

Сплайны соответствующей степени строятся аналогично. Наберем с клавиатуры

KLS := lspline(X ,Y ) KPS := pspline(X ,Y ) KCS := cspline(X ,Y )

y1(x):= interp(KLS, X ,Y , x) y2(x):= interp(KPS, X ,Y , x) y3(x):= interp(KCS, X ,Y , x).

На всем отрезке x [0.58,0.64] графики имеют вид

На концах отрезка вблизи точек x = 0.58 и x = 0.64 все четыре графика построим от-

дельно:
x := 0.56,0.5605…0.62	xx := 0.62,0.6205…0.64

Явно видимое различие трех графиков проявляется лишь при экстраполяции; внутри отрезка x [0.58,0.64] все три графика почти сливаются в один.

Из всех четырех графиков видно, что функция lspline(X,Y) строит вектор коэффициентов линейного сплайна дефекта нуль, то есть во всех внутренних узлах сетки состыкованы сама функция y = y(x) и ее первая производная. Линейный сплайн дефекта один – это ин-

терполяция ломаными, то есть функция linterp(X,Y,x), график которой приведен на самом первом рисунке.

К сожалению, в доступной автору документации пакета Mathcad нет упоминания, какого типа краевые условия употреблялись для каждого вида сплайнов. Скорее всего, это были так называемые естественные условия, то есть условия равенства нулю старшей используемой производной на концах рассматриваемого отрезка.

Составим программу вычисления коэффициентов кубического сплайна с краевыми условиями четвертого типа – «отсутствие узла». Эти коэффициенты вычисляются решением системы линейных уравнений (3.10.4) методом прогонки. Система (3.10.4) отличается от классической трехдиагональной системы уравнений, используемой в методе прогонки, тем,

что в первом и последнем уравнении этой системы имеется по одному лишнему (ненулевому) коэффициенту. Поэтому формулы метода прогонки, приведенные в подразд. 5.4, необходимо несколько модернизировать. Модернизация касается способа вычисления x1 и x2 в

прямом ходе и x1 в обратном; все остальные переменные определяются стандартно по фор-

мулам (5.4.1) – (5.4.6).

Progonka1–подпрограмма, реализующая модернизированный метод прогонки. Параметры: вектор b содержит коэффициенты при неизвестных левой части системы уравнений, стоящие на главной диагонали, вектор a -под главной диагональю, вектор c -над главной

диагональю. Коэффициент p1 первого уравнения системы b1 x1 + c1 x2	+ p1 x3	= d1 находится в
cn , коэффициент p2 последнего уравнения p2 xn−2 + an xn−1 + bn xn = dn	в a1 .	Вектор d содер-

жит свободные члены (правые части системы уравнений). Подпрограмма выдает вектор x — вектор решений системы. Для правильной работы подпрограммы необходимо задание

ORIGIN :=1.

Следующая подпрограмма intspline вычисляет значение кубического сплайна с краевыми условиями четвертого типа для данного значения аргумента x . X и Y — векторы исходных данных, n — число точек сетки. Структура подпрограммы intspline совершенно прозрачна. Сначала вычисляются все коэффициенты диагональной системы (3.10.4), затем эта система решается модифицированным методом прогонки. Далее по значению аргумента x

выбирается нужный кубический многочлен третьей степени P3 (x) и находится его значение

при данном x по формуле (3.10.1) кубического интерполяционного многочлена Эрмита.

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

KCS := cspline (X ,Y ) y3(x1):= interp(KCS , X ,Y , x1)

i := 1…100 x2i := X 1 + (X 7 − X 1 ) (i −1) 100

y1i := intspline (X, Y,7, x2i )

Видно, что оба графика на отрезке x [0.59,0.63] абсолютно совпадают.

В этом можно убедиться и сравнить полученные интерполированные значения f (x)

по двум разным программам:

x5 := 0.6 y5 := interp(KCS, X ,Y , x5) y5 = 0.536 y55 := intspline(X ,Y ,7, x5) y55 = 0.536

x6 := 0.63 y6 := interp(KCS, X ,Y , x6) y6 = 0.560 y56 := intspline(X ,Y ,7, x6) y56 = 0.560.

Задание № 1. Построить для таблично заданной функции y = f (x) кубический сплайн с граничными условиями любого из четырех типов на концах отрезка:

				Номера вариантов
	1		2		3		4		5
xi	yi	xi	yi	xi	yi	xi	yi	xi	yi
0.698	2.2234	0.100	1.1213	0.235	1.2080	0.095	1.0913	0.103	2.0128
0.706	2.2438	0.108	1.1316	0.240	1.2126	0.102	1.2349	0.108	2.0334
0.714	2.2645	0.119	1.1459	0.250	1.2217	0.104	1.2799	0.115	2.0607
0.727	2.2984	0.127	1.1565	0.255	1.2263	0.107	1.3514	0.120	2.0792
0.736	2.3222	0.135	1.1671	0.265	1.2355	0.110	1.4282	0.128	2.1072
0.747	2.3516	0.146	1.1819	0.280	1.2493	0.112	1.4826	0.136	2.1335
0.760	2.3869	0.157	1.1969	0.295	1.2633	0.116	1.6003	0.141	2.1492
0.769	2.4116	0.169	1.2134	0.300	1.2680	0.120	1.7321	0.150	2.1761
0.782	2.4478	0.175	1.2196	0.305	1.2726	0.125	1.8500	0.157	2.1805
				Номера вариантов
	6		7		8		9		10
xi	yi	xi	yi	xi	yi	xi	yi	xi	yi
0.296	3.2558	0.050	0.2079	0.902	1.2351	0.100	1.8378	0.400	1.6683
0.303	3.1764	0.052	0.2081	0.909	1.2369	0.104	1.8369	0.405	1.6664
0.310	3.1218	0.060	0.2090	0.919	1.2394	0.118	1.8335	0.410	1.6645
0.317	3.0482	0.065	0.2095	0.940	1.2448	0.139	1.8286	0.420	1.6607
0.323	2.9876	0.069	0.2099	0.944	1.2458	0.145	1.8272	0.429	1.6573
0.330	2.9195	0.075	0.2105	0.955	1.2486	0.158	1.8242	0.440	1.6532
0.339	2.8360	0.085	0.2116	0.965	1.2511	0.167	1.8221	0.449	1.6499
0.345	2.7771	0.090	0.2121	0.975	1.2537	0.185	1.8179	0.455	1.6476
0.352	2.6114	0.096	0.2127	1.010	1.2628	0.200	1.8145	0.465	1.6439
				Номера вариантов
	11		12		13		14		15
xi	yi	xi	yi	xi	yi	xi	yi	xi	yi
1.030	2.8011	0.200	0.0400	0.010	1.0101	0.100	0.8100	1.000	0.0000
1.080	2.9447	0.300	0.0899	0.105	1.1105	0.155	0.7140	1.510	0.2729
1.160	3.1899	0.350	0.1222	0.156	1.1681	0.220	0.6084	2.100	0.3533
1.230	3.4212	0.379	0.1431	0.200	1.2198	0.280	0.5184	2.750	0.3679
1.260	3.5254	0.415	0.1714	0.215	1.2378	0.375	0.3906	3.420	0.3595
1.330	3.7810	0.500	0.2474	0.289	1.3298	0.445	0.3080	3.915	0.3486
1.390	4.0149	0.596	0.3478	0.316	1.3645	0.510	0.2401	4.350	0.3380
1.450	4.1713	0.615	0.3693	0.390	1.4626	0.625	0.1406	4.800	0.3268
1.500	4.2390	0.700	0.4706	0.500	1.6152	0.770	0.0529	5.200	0.3171
				Номера вариантов
	16		17		18		19		20
xi	yi	xi	yi	xi	yi	xi	yi	xi	yi
0.100	0.0010	-2.100	-0.674	-1.000	1.368	1.000	0.000	1.000	0.6663
0.210	0.0441	-1.815	0.914	-0.010	1.990	5.000	0.999	1.300	0.5706
0.295	0.0868	-1.600	0.574	0.750	3.117	15.150	0.411	1.650	0.5429
0.348	0.1205	-1.420	0.827	1.900	7.686	27.900	-0.186	1.915	0.5887
0.419	0.1738	-1.290	0.931	2.815	17.692	35.100	-0.405	2.280	0.7256
0.475	0.2219	-0.750	0.983	3.100	23.198	44.000	-0.599	2.500	0.8262
					95

Соседние файлы в предмете Вычислительная математика

Источник

Макеты страниц

функция (рис. 11.8), являющаяся сплайном первой степени (линейным сплайном) с дефектом, равным единице. Действительно, на отрезке сама функция (нулевая производная) непрерывна. В то же время на каждом частичном отрезке совпадает с некоторым многочленом первой степени.

Рис. 11.8

Рис. 11.9

Наиболее широкое распространение на практике получили сплайны третьей степени (кубические сплайны) с дефектом, равным 1 или 2. Такие сплайны на каждом из частичных отрезков совпадают с кубическим многочленом:

и имеют на отрезке по крайней мере одну непрерывную производную

Термин «сплайн» происходит от английского слова (гибкая линейка, стержень) — названия приспособления, использовавшегося чертежниками для проведения гладких кривых через заданные точки. Бели гибкую стальную линейку поставить на ребро и, изогнув, зафиксировать ее положение в узловых точках (рис. 11.9), то получится механический аналог кубического сплайна. В самом деле, из курса сопротивления материалов известно, что уравнение свободного равновесия профиля линейки таково: Следовательно, в промежутке между двумя соседними узлами представляет собой многочлен третьей степени. В то же время отсутствие у линейки изломов свидетельствует о непрерывности касательной к графику функции и кривизны, т. е. производных

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

Пусть функция задана таблицей своих значений Сплайн называется интерполяционным, если для всех Значение называется наклоном сплайна в точке

Заметим, что на отрезке интерполяционный кубический сплайн однозначно определяется заданием значений самом деле, из равенства (11.31) вытекает следующая формула:

Здесь

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

3. Локальный сплайн.

Если в точках х, известны значения производной то естественно положить для всех Тогда на каждом частичном отрезке в соответствии с формулой (11.64) сплайн однозначно определяется значениями (поэтому его и называют локалъныл сплайном). Заметим, что он совпадает с кубическим интерполяционным многочленом Эрмита (11.31) для отрезка

Из неравенства (11.33) получается следующая оценка погрешности интерполяции локальным кубическим сплайном:

где Ащах максимальная из длин частичных отрезков.

Заметим, что для построенного указанным образом сплайна можно гарантировать непрерывность на отрезке только функции и ее первой производной 53, т.е. его дефект равен 2.

Существуют и другие способы выбора коэффициентов а, приводящие к локальным сплайнам (кубический многочлен Бесселя, метод Акимы и др. [16]).

4. Глобальные способы построения кубических сплайнов.

Для того чтобы сплайн имел непрерывную на отрезке вторую производную необходимо выбирать наклоны а, так, чтобы в точках х, «стыка» многочленов совпадали значения их вторых производных:

Пользуясь формулой (11.64), найдем значение

Из подобной формулы, записанной

Таким образом, равенства (11.66) приводят к следующей системе уравнений относительно коэффициентов

Заметим, что эта система уравнений недоопределена, так как число уравнений системы (равное меньше числа неизвестных (равного Выбор двух оставшихся уравений обычно связывают с некоторыми дополнительными условиями, накладываемыми на сплайн в граничных точках (граничными условиями). Укажем на некоторые из наиболее известных граничных условий.

1°. Если в граничных точках известны значения первой производной то естественно положить

Дополняя систему (11.69) уравнениями (11.70), приходим к системе уравнений с трехдиагональной матрицей, которая легко решается методом прогонки (см. гл. 5). Полученный таким образом сплайн называется фундаментальным кубическим сплайном.

2°. Если в граничных точках известны значения второй производной то можно наложить на сплайн граничные условия что приводит к следующим уравнениям:

(достаточно в равенстве (11.68) взять а в равенстве

3°. Полагая в уравнениях от того, выполнены ли эти условия для интерполируемой функции), придем к системе уравнений, определяющих так называемый естественный кубический сплайн.

4°. Часто нет, никакой дополнительной информации о значениях производных на концах отрезка. Один из применяемых в этой ситуациии подходов состоит в использовании условия «отсутствия узла». Выбор наклонов производят таким образом, чтобы для получаемого сплайна выполнялись условия Для этого достаточно потребовать совпадения в точках соответствующих третьих производных:

Эквивалентные алгебраические уравнения выглядят так:

Та же аппроксимирующая функция может быть получена несколько иначе. Уменьшим число частичных отрезков, объединив попарно отрезки Это отвечает разбиению отрезка точками где для и построению соответствующего интерполяционного сплайна Условия «отсутствия узла» эквидалентны требованию совпадения значений сплайна в точках со значениями

5°. Если периодическая функция с периодом, равным о то систему (11-69) следует дополнить условиями

Существуют и другие подходы к заданию граничных условий (подробнее об этом см. [16]).

Пример 11.13. Для функции, заданной табл. 11.12, построй (естественный) кубический сплайн. В этом случае система уравнений для наклонов в точках записывается следующим образом:

Решал ее, получаем значения . Теперь на каждом частичном отрезке значения сплайна можно вычислить по формуле (11.64). Соответствующий график приведен на рис 11.10 (ср. с рис. 11.6 и 11.7).

Пример 11.14. Интерполируем функцию, заданную табл. 11.12, кубическим сплайном, используя условие «отсутствия узла». В этом случае уравнения останутся прежними, а уравнения (11 75) и (11.79) заменяются следующими:

Рис. 11.10

Решая систему получаем значения График соответствующего сплайна мало отличается от графика, изображенного на рис. 11.10.

5. Погрешность приближения кубическими сплайнами.

Теорема 11.9. Пусть функция имеет на отрезке непрерывную производную четвертою порядка и Тогда для интерполяционною кубическою сплайна удовлетворяющею граничным условиям типов 1°, 2°, 4° иди 5° (последнее — для случая периодической функции), справедлива следующая оценка погрешности:

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

Теорема 11.10. При выполнении условий теоремы 11.9 для указанных в ней сплайнов справедливы неравенства

Замечание. Благодаря большей простоте записи и благозвучному названию естественные сплайны получили значительное распространение. Однако искусственное наложение условий при интерполяции функций, которые этим условиям не удовлетворяют, приводит к значительной потере точности. Вместо четвертого порядка точности (как локальный кубический сплайн или кубические сплайны с граничными условиями типов 1°, 2°, 4°, 5°) естественный сплайн обладает лишь вторым порядком точности. Если использование естественного сплайна не вызвано какими-либо специальными причинами, то следует, по-видимому, отказаться от него в пользу кубического сплайна с граничным условием типа 4°.

Оглавление

ПРЕДИСЛОВИЕ
Глава 1. МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ И РЕШЕНИЕ ИНЖЕНЕРНЫХ ЗАДАЧ С ПРИМЕНЕНИЕМ ЭВМ
§ 1.2. Основные этапы решения инженерной задачи с применением ЭВМ
§ 1.3. Вычислительный эксперимент
§ 1.4. Дополнительные замечания
Глава 2. ВВЕДЕНИЕ В ЭЛЕМЕНТАРНУЮ ТЕОРИЮ ПОГРЕШНОСТЕЙ
§ 2.1. Источники и классификация погрешностей результата численного решения задачи
§ 2.2. Приближенные числа. Абсолютная и относительная погрешности
2. Правила записи приближенных чисел.
3. Округление.
§ 2.3. Погрешности арифметических операций над приближенными числами
§ 2.4. Погрешность функции
§ 2.5. Особенности машинной арифметики
2. Представление целых чисел.
3. Представление вещественных чисел.
4. Арифметические операции над числами с плавающей точкой.
5. Удвоенная точность.
6. Вычисление машинного эпсилон.
§ 2.6. Дополнительные замечания
Глава 3. ВЫЧИСЛИТЕЛЬНЫЕ ЗАДАЧИ, МЕТОДЫ И АЛГОРИТМЫ. ОСНОВНЫЕ ПОНЯТИЯ
§ 3.2. Обусловленность вычислительной задачи
2. Примеры плохо обусловленных задач.
3. Обусловленность задачи вычисления значения функции одной переменной.
4. Обусловленность задачи вычисления интеграла …
5. Обусловленность задачи вычисления суммы ряда.
§ 3.3. Вычислительные методы
§ 3.4. Корректность вычислительных алгоритмов
§ 3.5. Чувствительность вычислительных алгоритмов к ошибкам округления
§ 3.6. Различные подходы к анализу ошибок
§ 3.7. Требования, предъявляемые к вычислительным алгоритмам
§ 3.8. Дополнительные замечания
Глава 4. МЕТОДЫ ОТЫСКАНИЯ РЕШЕНИЙ НЕЛИНЕЙНЫХ УРАВНЕНИЙ
§ 4.2. Обусловленность задачи вычисления корня
§ 4.3. Метод бисекции
§ 4.4. Метод простой итерации
§ 4.5. Обусловленность метода простой итерации
§ 4.6. Метод Ньютона
§ 4.7. Модификации метода Ньютона
§ 4.8. Дополнительные замечания
Глава 5. ПРЯМЫЕ МЕТОДЫ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ
§ 5.2. Нормы вектора и матрицы
§ 5.3. Типы используемых матриц
§ 5.4. Обусловленность задачи решения системы линейных алгебраических уравнений
§ 5.5 Метод Гаусса
§ 5.6. Метод Гаусса и решение систем уравнений с несколькими правыми частями, обращение матриц, вычисление определителей
§ 5.7. Метод Гаусса и разложение матрицы на множители. LU-разложение
§ 5.8. Метод Холецкого (метод квадратных корней)
§ 5.9. Метод прогонки
§ 5.10. QR-разложение матрицы. Методы вращений и отражений
§ 5.11. Итерационное уточнение
§ 5.12. Дополнительные замечания
Глава 6. ИТЕРАЦИОННЫЕ МЕТОДЫ РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ
§ 6.1. Метод простой итерации
§ 6.2. Метод Зейделя
§ 6.3. Метод релаксации
§ 6.4. Дополнительные замечания
Глава 7. МЕТОДЫ ОТЫСКАНИЯ РЕШЕНИЙ СИСТЕМ НЕЛИНЕЙНЫХ УРАВНЕНИЙ
§ 7.2. Метод простой итерации
§ 7.3. Метод Ньютона для решения систем нелинейных уравнений
7.4. Модификации метода Ньютона
§ 7.5. О некоторых подходах к решению задач локализации и отыскания решений систем нелинейных уравнений
§ 7.6. Дополнительные замечания
Глава 8. МЕТОДЫ РЕШЕНИЯ ПРОБЛЕМЫ СОБСТВЕННЫХ ЗНАЧЕНИЙ
§ 8.2. Степенной метод
§ 8.3. Метод обратных итераций
§ 8.4. QR-алгоритм
§ 8.5. Дополнительные замечания
Глава 9. МЕТОДЫ ОДНОМЕРНОЙ МИНИМИЗАЦИИ
§ 9.2. Обусловленность задачи минимизации
§ 9.3. Методы прямого поиска. Оптимальный пассивный поиск. Метод деления отрезка пополам. Методы Фибоначчи и золотого сечения
§ 9.4. Метод Ньютона и другие методы минимизация гладких функций
§ 9.5. Дополнительные замечания
Глава 10. МЕТОДЫ МНОГОМЕРНОЙ МИНИМИЗАЦИИ
§ 10.1. Задача безусловной минимизации функции многих переменных
§ 10.2. Понятие о методах спуска. Покоординатный спуск
§ 10.3. Градиентный метод
§ 10.4. Метод Ньютона
§ 10.5. Метод сопряженных градиентов
§ 10.6. Метода минимизации без вычисления производных
§ 10.7. Дополнительные замечания
Глава 11. ПРИБЛИЖЕНИЕ ФУНКЦИЙ И СМЕЖНЫЕ ВОПРОСЫ
§ 11.2. Интерполяция обобщенными многочленами
§ 11.3. Полиномиальная интерполяция. Многочлен Лагранжа
§ 11.4. Погрешность интерполяции
§ 11.5. Интерполяция с кратными узлами
§ 11.6. Минимизация оценки погрешности интерполяции. Многочлены Чебышева
§ 11.7. Конечные разности
§ 11.8. Разделенные разности
§ 11.9. Интерполяционный многочлен Ньютона. Схема Эйткена
§ 11.10. Обсуждение глобальной полиномиальной интерполяции. Понятие о кусочно-полиномиальной интерполяции
§ 11.11. Интерполяция сплайнами
§ 11.12. Понятие о дискретном преобразовании Фурье и тригонометрической интерполяции
§ 11.13. Метод наименьших квадратов
§ 11.14. Равномерное приближение функций
§ 11.15. Дробно-рациональные аппроксимации и вычисление элементарных функций
§ 11.16. Дополнительные замечания
Глава 12. ЧИСЛЕННОЕ ДИФФЕРЕНЦИРОВАНИЕ
§ 12.1. Простейшие формулы численного дифференцирования
§ 12.2. О выводе формул численного дифференцирования
§ 12.3. Обусловленность формул численного дифференцирования
§ 12.4. Дополнительные замечания
Глава 13. ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ
13.2. Квадратурные формулы интерполяционного типа
§ 13.3. Квадратурные формулы Гаусса
§ 13.4. Апостериорные оценки погрешности. Понятие об адаптивных процедурах численного интегрирования
§ 13.5. Вычисление интегралов в нерегулярных случаях
§ 13.6. Дополнительные замечания
Глава 14. ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ ЗАДАЧИ КОШИ ДЛЯ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ
§ 14.1. Задача Коши для дифференциального уравнения первого порядка
§ 14.2. Численные методы решения задачи Коши. Основные понятия и определения
§ 14.3. Использование формулы Тейлора
§ 14.4. Метод Эйлера
§ 14.5. Модификации метода Эйлера второго порядка точности
§ 14.6. Методы Рунге-Кутты
§ 14.7. Линейные многошаговые методы. Методы Адамса
§ 14.8. Устойчивость численных методов решения задачи Коши
§ 14.9. Неявный метод Эйлера
§ 14.10. Решение задачи Коши для систем обыкновенных дифференциальных уравнений и дифференциальных уравнений m-го порядка
§ 14.11. Жесткие задачи
§ 14.12. Дополнительные замечания
Глава 15. РЕШЕНИЕ ДВУХТОЧЕЧНЫХ КРАЕВЫХ ЗАДАЧ
§ 15.1. Краевые задачи для одномерного стационарного уравнения теплопроводности
§ 15.2. Метод конечных разностей: основные понятия
§ 15.3. Метод конечных разностей: аппроксимации специального вида
§ 15.4. Понятие о проекционных и проекционно-разностных методах. Методы Ритца и Гадеркина. Метод конечных элементов
§ 15.5. Метод пристрелки
§ 15.6. Дополнительные замечания

Источник

This example shows how to use the csapi and csape commands from Curve Fitting Toolbox™ to construct cubic spline interpolants.

The CSAPI Command

The command

values = csapi(x,y,xx)

returns the values at xx of the cubic spline interpolant to the given data (x,y), using the not-a-knot end condition. This interpolant is a piecewise cubic function, with break sequence x, whose cubic pieces join together to form a function with two continuous derivatives. The «not-a-knot» end condition means that, at the first and last interior break, even the third derivative is continuous (up to round-off error).

Specifying only two data points results in a straight line interpolant.

x = [0 1];
y = [2 0];
xx = linspace(0,6,121);
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Interpolant to Two Points')

Figure contains an axes object. The axes object with title Interpolant to Two Points contains 2 objects of type line.

Specifying three data points gives a parabola.

x = [2 3 5];
y = [1 0 4];
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Interpolant to Three Points')

Figure contains an axes object. The axes object with title Interpolant to Three Points contains 2 objects of type line.

More generally, four or more data points give a cubic spline.

x = [1 1.5 2 4.1 5];
y = [1 -1 1 -1 1];
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Cubic Spline Interpolant to Five Points')

Figure contains an axes object. The axes object with title Cubic Spline Interpolant to Five Points contains 2 objects of type line.

How to Check the Output of CSAPI

These look like nice interpolants, but how do we check that csapi performs as advertised?

We already saw that csapi interpolates, because we plotted the data points and the interpolant went right through those points. But to be sure that we get a cubic spline, it is best to start with data from a cubic spline of the expected sort and check whether csapi reproduces that cubic spline, i.e., gives back that cubic spline from which the data were taken.

Example: The Truncated Power Function

One simple example of a cubic spline function to check against is the truncated third power, i.e., the function

f(x)=((x-xi)+)3,

where xi is one of the breaks and the «+» subscript indicates the truncation function, provided by the command subplus:

 SUBPLUS Positive part.
 
                                   x , if  x>=0
    y  = subplus(x) := (x)_{+}  =               ,
                                   0 , if  x<=0
 
   returns the positive part of X. Used for computing truncated powers.

The truncated 3rd power is plotted below for the particular choice xi = 2. As expected, it is zero to the left of 2, and rises like (x-2)^3 to the right of 2.

plot(xx, subplus(xx-2).^3,'y','LineWidth',3)
axis([0,6,-10,70])

Figure contains an axes object. The axes object contains an object of type line.

Now we interpolate this particular cubic spline at the data sites 0:6, and plot the interpolant on top of the spline, in black.

x = 0:6;
y = subplus(x-2).^3;
values = csapi(x,y,xx);
hold on
plot(xx,values,'k',x,y,'ro')
hold off
title('Interpolant to ((x-2)_+)^3')

Figure contains an axes object. The axes object with title I n t e r p o l a n t blank t o blank ( ( x - 2 ) indexOf + baseline ) Cubed baseline contains 3 objects of type line.

The Interpolation Error

When comparing two functions, it is usually much more informative to plot their difference.

plot(xx, values - subplus(xx-2).^3)
title('Error in Cubic Spline Interpolation to ((x-2)_+)^3')

Figure contains an axes object. The axes object with title E r r o r blank i n blank C u b i c blank S p l i n e blank I n t e r p o l a t i o n blank t o blank ( ( x - 2 ) indexOf + baseline ) Cubed baseline contains an object of type line.

To put the size of their difference into context, you can also compute the maximum data value. This shows the error to be no worse than the inevitable round-off error.

A Truncated Power That Cannot be Reproduced

As a further test, we interpolate a truncated power whose csapi-produced interpolant at the sites 0:6 cannot coincide with it. For example, the first interior break of the interpolating spline is not really a knot since csapi uses the «not-a-knot» condition, hence the interpolant has three continuous derivatives at that site. This implies that we should not be able to reproduce the truncated 3rd power centered at that site since its third derivative is discontinuous across that site.

values = csapi(x,subplus(x-1).^3,xx);
plot(xx, values - subplus(xx-1).^3)
title('Error in Not-a-Knot Interpolant to ((x-1)_+)^3')

Figure contains an axes object. The axes object with title E r r o r blank i n blank N o t - a - K n o t blank I n t e r p o l a n t blank t o blank ( ( x - 1 ) indexOf + baseline ) Cubed baseline contains an object of type line.

Since 1 is a first interior knot, it is not active for this interpolant.

The difference is as large as .18, but decays rapidly as we move away from 1. This illustrates that cubic spline interpolation is essentially local.

Using the ppform Instead of Values

It is possible to retain the interpolating cubic spline in a form suitable for subsequent evaluation, or for calculating its derivatives, or for other manipulations. This is done by calling csapi in the form

pp = csapi(x,y)

which returns the ppform of the interpolant. You can evaluate this form at some new points xx by the command

values = fnval(pp,xx)

You can differentiate the interpolant by the command

dpp = fnder(pp)

or integrate it by the command

ipp = fnint(pp)

which return the ppform of the derivative or the integral, respectively.

Example: Differentiating and Integrating the Interpolant

To show differentiation of an interpolant, we plot the derivative of this truncated power

f2′(x)=3((x-2)+)2,

(again in yellow) and then, on top of it, the derivative of our interpolant to the original truncated third power function (again in black).

plot(xx,3*subplus(xx-2).^2,'y','LineWidth',3)
pp = csapi(x,subplus(x-2).^3);
dpp = fnder(pp);
hold on
plot(xx,fnval(dpp,xx),'k')
hold off
title('Derivative of Interpolant to ((x-2)_+)^3')

Figure contains an axes object. The axes object with title D e r i v a t i v e blank o f blank I n t e r p o l a n t blank t o blank ( ( x - 2 ) indexOf + baseline ) Cubed baseline contains 2 objects of type line.

Again, the more informative comparison is to plot their difference, and as before this is no bigger than the round-off error.

plot(xx, fnval(dpp,xx) - 3*subplus(xx-2).^2)
title('Error in Derivative of interpolant to ((x-2)_+)^3')

Figure contains an axes object. The axes object with title E r r o r blank i n blank D e r i v a t i v e blank o f blank i n t e r p o l a n t blank t o blank ( ( x - 2 ) indexOf + baseline ) Cubed baseline contains an object of type line.

The second derivative of the truncated power is

f2′′(x)=6(x-2)+

A plot of the difference between this function and the second derivative of the interpolant to the original function shows that there are now jumps, but they are still within the round-off error.

ddpp = fnder(dpp);
plot(xx, fnval(ddpp,xx) - 6*subplus(xx-2))
title('Error in Second Derivative of Interpolant to ((x-2)_+)^3')

Figure contains an axes object. The axes object with title E r r o r blank i n blank S e c o n d blank D e r i v a t i v e blank o f blank I n t e r p o l a n t blank t o blank ( ( x - 2 ) indexOf + baseline ) Cubed baseline contains an object of type line.

The integral of the truncated power is

F2(x)=((x-2)+)4/4.

A plot of the difference between this function and the integral of the interpolant to the original function again shows that the errors are within the round-off error.

ipp = fnint(pp);
plot(xx, fnval(ipp,xx) - subplus(xx-2).^4/4)
title('Error in Integral of Interpolant to ((x-2)_+)^3')

Figure contains an axes object. The axes object with title E r r o r blank i n blank I n t e g r a l blank o f blank I n t e r p o l a n t blank t o blank ( ( x - 2 ) indexOf + baseline ) Cubed baseline contains an object of type line.

The CSAPE Command

Like csapi, the csape command provides a cubic spline interpolant to given data. However, it permits various additional end conditions. Its simplest version,

pp = csape(x,y)

uses the Lagrange end condition, which is a common alternative to the not-a-knot condition used by csapi. csape does not directly return values of the interpolant, but only its ppform.

For example, consider again interpolation to the function

f1(x)=((x-1)+)3,

which csapi fails to reproduce well. We plot the error of the not-a-knot interpolant returned by csapi (in black), along with the error of the interpolant obtained from csape (in red).

exact = subplus(xx-1).^3;
plot(xx, fnval(csapi(x,subplus(x-1).^3),xx) - exact,'k')
hold on
plot(xx, fnval(csape(x,subplus(x-1).^3),xx) - exact,'r')
title('Error in Not-a-Knot vs. Lagrange End Conditions')
legend({'Not-a-Knot' 'Lagrange'});
hold off

Figure contains an axes object. The axes object with title Error in Not-a-Knot vs. Lagrange End Conditions contains 2 objects of type line. These objects represent Not-a-Knot, Lagrange.

There is not much difference between the two interpolants in this case.

Other End Conditions: The ‘Natural’ Spline Interpolant

The csape command also provides ways to specify several other types of end conditions for an interpolating cubic spline. For example, the command

pp = csape(x,y,'variational')

uses the so-called ‘natural’ end conditions. This means that the second derivative is zero at the two extreme breaks.

This step shows how to apply ‘natural’ cubic spline interpolation to the function

f2(x)=((x-2)+)3,

and plot the error. The code below computes the ‘natural’ spline interpolant with an alternative argument syntax that is equivalent to the 'variational' argument: using 'second' specifies that csape should set the second derivative at the extreme data sites to the default value of 0.

pp = csape(x,subplus(x-2).^3,'second');
plot(xx, fnval(pp,xx) - subplus(xx-2).^3)
title('Error in ''Natural'' Spline Interpolation to ((x-2)_+)^3')

Figure contains an axes object. The axes object with title E r r o r blank i n blank ' N a t u r a l ' blank S p l i n e blank I n t e r p o l a t i o n blank t o blank ( ( x - 2 ) indexOf + baseline ) Cubed baseline contains an object of type line.

Note the large error near the right end. This is due to the fact that the ‘natural’ end conditions implicitly insist on having a zero second derivative there.

Other End Conditions: Prescribing Second Derivatives

We can also explicitly use the correct second derivatives to get a small error. First, we compute the correct second derivative values of the truncated power at the endpoints.

endcond = 6*subplus(x([1 end])-2);

Then we create the interpolant, specifying that second derivatives at the endpoints are to be matched to the second derivative values we just computed. We do this by providing endcond(1) for the left endpoint condition, and endcond(2) for the right, along with the data values.

pp = csape(x,[endcond(1) subplus(x-2).^3 endcond(2)], 'second');
plot(xx, fnval(pp,xx) - subplus(xx-2).^3,'r')
title(['Error in Spline Interpolation to ((x-1)_+)^3'; ...
       '  When Matching the 2nd Derivative at Ends  '])

Figure contains an axes object. The axes object with title E r r o r blank i n blank S p l i n e blank I n t e r p o l a t i o n blank t o blank ( ( x - 1 ) indexOf + baseline ) Cubed baseline blank blank blank W h e n blank M a t c h i n g blank t h e blank 2 n d blank D e r i v a t i v e blank a t blank E n d s blank blank contains an object of type line.

Other End Conditions: Prescribing Slopes

csape also permits specification of endpoint slopes. This is the clamped (or, complete) cubic spline interpolant. The statement

pp = csape(x,[sl,y,sr],'clamped')

creates the cubic spline interpolant to the data (x, y) that also has slope sl at the leftmost data site and slope sr at the rightmost data site.

Other End Conditions: Mixed End Conditions

It is even possible to mix these conditions. For example, our much-exercised truncated power function

f1(x)=((x-1)+)3

has slope 0 at x=0 and second derivative 30 at x=6 (the last data site).

Therefore, by matching the slope at the left end and the curvature at the right, we expect no error in the resulting interpolant.

pp = csape(x, [0 subplus(x-1).^3 30], [1 2]);
plot(xx, fnval(pp,xx) - subplus(xx-1).^3)
title(['Error in Spline Interpolation to ((x-1)_+)^3'; ...
       '        with Mixed End Conditions.          '])

Other End Conditions: Periodic Conditions

It is also possible to prescribe periodic end conditions. For example, the sine function is 2*pi-periodic and has the values [0 -1 0 1 0] at the sites (pi/2)*(-2:2). The difference, between the sine function and its periodic cubic spline interpolant at these sites, is only 2 percent. Not bad.

x = (pi/2)*(-2:2);
y = [0 -1 0 1 0];
pp = csape(x,y, 'periodic' );
xx = linspace(-pi,pi,201);
plot(xx, sin(xx) - fnval(pp,xx), 'x')
title('Error in Periodic Cubic Spline Interpolation to sin(x)')

Figure contains an axes object. The axes object with title Error in Periodic Cubic Spline Interpolation to sin(x) contains an object of type line.

End Conditions Not Explicitly Covered by CSAPI or CSAPE

Any end condition not covered explicitly by csapi or csape can be handled by constructing the interpolant with the csape default side conditions, and then adding to it an appropriate scalar multiple of an interpolant to zero values and some side conditions. If there are two `nonstandard’ side conditions to be satisfied, you may have to solve a 2-by-2 linear system first.

For example, suppose that you want to compute the cubic spline interpolant s to the data

x = 0:.25:3;
q = @(x) x.*(-1 + x.*(-1+x.*x/5));
y = q(x);

and enforce the condition

lambda(s) := a * (Ds)(e) + b * (D^2 s)(e) = c

on the first and second derivatives of s at the point e.

The data were generated from a quartic polynomial that happens to satisfy this side condition with specific parameters

e = x(1);
a = 2; b = -3; c = 4;

To construct the interpolant that satisfies this specific condition, we first construct the interpolant with the default end conditions

and the first derivative of its first polynomial piece.

dp1 = fnder(fnbrk(pp1,1));

In addition, we construct the cubic spline interpolant to zero data values, specifying that it have a slope of 1 at e,

pp0 = csape(x,[1,zeros(size(y)),0], [1,0]);

as well as constructing the first derivative of its first polynomial piece.

dp0 = fnder(fnbrk(pp0,1));

Then we compute lambda for both pp1 and pp0,

lam1 = a*fnval(dp1,e) + b*fnval(fnder(dp1),e);
lam0 = a*fnval(dp0,e) + b*fnval(fnder(dp0),e);

and construct the correct linear combination of pp1 and pp0 to get a cubic spline

s := pp1 + ((c - lambda(pp1))/lambda(pp0)) * pp0

that does satisfy the desired condition, as well as the default end condition at the right endpoint. We form this linear combination with the help of fncmb.

s = fncmb(pp0,(c-lam1)/lam0,pp1);

A plot of the interpolation error shows that s fits the quartic polynomial slightly better near e than the interpolant pp1 with the default conditions does.

xx = (-.3):.05:.7; yy = q(xx);
plot(xx, fnval(pp1,xx) - yy, 'x')
hold on
plot(xx, fnval(s,xx) - yy, 'o')
hold off
legend({'Default conditions' 'Nonstandard conditions'},'location','SE')

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent Default conditions, Nonstandard conditions.

If we want to enforce the condition

mu(s) := (D^3 s)(3) = 14.6

on the third derivative of the interpolant (the quartic satisfies this condition), then we construct an additional cubic spline interpolating to zero values, and with zero first derivative at the left endpoint, hence certain to be independent from pp0.

pp2 = csape(x,[0,zeros(size(y)),1],[0,1]);

Then we find the coefficients d0 and d2 in the linear combination

s := pp1 + d0*pp0 + d2*pp2

that solves the linear system

lambda(s) = c

mu(s) = 14.6

Note that both pp0 and pp2 vanish at all interpolation sites, hence s will match the given data for any choice of d0 and d2.

For amusement, we use the MATLAB® encoding facility to write a loop to compute lambda(pp_j) and mu(pp_j), for j=0:2.

dd = zeros(2,3);
for j=0:2
   J = num2str(j);
   eval(['dpp',J,'=fnder(pp',J,');']);
   eval(['ddpp',J,'=fnder(dpp',J,');']);
   eval(['dd(1,1+',J,')=a*fnval(dpp',J,',e)+b*fnval(ddpp',J,',e);']);
   eval(['dd(2,1+',J,')=fnval(fnder(ddpp',J,'),3);']);
end

Given the values of lambda and mu for pp0, pp1, and pp2, we then solve for the coefficients that define the correct linear combination.

d = dd(:,[1,3])([c;14.6]-dd(:,2));
s = fncmb(fncmb(pp0,d(1),pp2,d(2)),pp1);

xxx = 0:.05:3;
yyy = q(xxx);
plot(xxx, yyy - fnval(s,xxx),'x')
title('Error in Spline Interpolant to y = x*(-1 + x*(-1+x*x/5))')

Figure contains an axes object. The axes object with title Error in Spline Interpolant to y = x*(-1 + x*(-1+x*x/5)) contains an object of type line.

For reassurance, we compare this error with the one obtained in complete cubic spline interpolation to this function:

hold on
plot(xxx, yyy - fnval(csape(x,[-1,y,-7+(4/5)*27],'clamped'),xxx),'o')
hold off
legend({'Nonstandard conditions' 'Endslope conditions'})

Figure contains an axes object. The axes object with title Error in Spline Interpolant to y = x*(-1 + x*(-1+x*x/5)) contains 2 objects of type line. These objects represent Nonstandard conditions, Endslope conditions.

The errors differ (and not by much) only near the end points, testifying to the fact that both pp0 and pp2 are sizable only near their respective end points.

As a final check, we verify that s satisfies the desired third derivative condition at 3.

Источник

This example shows how to use the csapi and csape commands from Curve Fitting Toolbox™ to construct cubic spline interpolants.

The CSAPI Command

The command

values = csapi(x,y,xx)

Specifying only two data points results in a straight line interpolant.

x = [0 1];
y = [2 0];
xx = linspace(0,6,121);
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Interpolant to Two Points')

Figure contains an axes object. The axes object with title Interpolant to Two Points contains 2 objects of type line.

Specifying three data points gives a parabola.

x = [2 3 5];
y = [1 0 4];
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Interpolant to Three Points')

Figure contains an axes object. The axes object with title Interpolant to Three Points contains 2 objects of type line.

More generally, four or more data points give a cubic spline.

x = [1 1.5 2 4.1 5];
y = [1 -1 1 -1 1];
plot(xx,csapi(x,y,xx),'k-',x,y,'ro')
title('Cubic Spline Interpolant to Five Points')

Figure contains an axes object. The axes object with title Cubic Spline Interpolant to Five Points contains 2 objects of type line.

How to Check the Output of CSAPI

These look like nice interpolants, but how do we check that csapi performs as advertised?

Example: The Truncated Power Function

One simple example of a cubic spline function to check against is the truncated third power, i.e., the function

f(x)=((x-xi)+)3,

where xi is one of the breaks and the «+» subscript indicates the truncation function, provided by the command subplus:

 SUBPLUS Positive part.
 
                                   x , if  x>=0
    y  = subplus(x) := (x)_{+}  =               ,
                                   0 , if  x<=0
 
   returns the positive part of X. Used for computing truncated powers.

The truncated 3rd power is plotted below for the particular choice xi = 2. As expected, it is zero to the left of 2, and rises like (x-2)^3 to the right of 2.

plot(xx, subplus(xx-2).^3,'y','LineWidth',3)
axis([0,6,-10,70])

Figure contains an axes object. The axes object contains an object of type line.

Now we interpolate this particular cubic spline at the data sites 0:6, and plot the interpolant on top of the spline, in black.

x = 0:6;
y = subplus(x-2).^3;
values = csapi(x,y,xx);
hold on
plot(xx,values,'k',x,y,'ro')
hold off
title('Interpolant to ((x-2)_+)^3')

Figure contains an axes object. The axes object with title I n t e r p o l a n t blank t o blank ( ( x - 2 ) indexOf + baseline ) Cubed baseline contains 3 objects of type line.

The Interpolation Error

When comparing two functions, it is usually much more informative to plot their difference.

plot(xx, values - subplus(xx-2).^3)
title('Error in Cubic Spline Interpolation to ((x-2)_+)^3')

To put the size of their difference into context, you can also compute the maximum data value. This shows the error to be no worse than the inevitable round-off error.

A Truncated Power That Cannot be Reproduced

values = csapi(x,subplus(x-1).^3,xx);
plot(xx, values - subplus(xx-1).^3)
title('Error in Not-a-Knot Interpolant to ((x-1)_+)^3')

Since 1 is a first interior knot, it is not active for this interpolant.

The difference is as large as .18, but decays rapidly as we move away from 1. This illustrates that cubic spline interpolation is essentially local.

Using the ppform Instead of Values

pp = csapi(x,y)

which returns the ppform of the interpolant. You can evaluate this form at some new points xx by the command

values = fnval(pp,xx)

You can differentiate the interpolant by the command

dpp = fnder(pp)

or integrate it by the command

ipp = fnint(pp)

which return the ppform of the derivative or the integral, respectively.

Example: Differentiating and Integrating the Interpolant

To show differentiation of an interpolant, we plot the derivative of this truncated power

f2′(x)=3((x-2)+)2,

(again in yellow) and then, on top of it, the derivative of our interpolant to the original truncated third power function (again in black).

plot(xx,3*subplus(xx-2).^2,'y','LineWidth',3)
pp = csapi(x,subplus(x-2).^3);
dpp = fnder(pp);
hold on
plot(xx,fnval(dpp,xx),'k')
hold off
title('Derivative of Interpolant to ((x-2)_+)^3')

Again, the more informative comparison is to plot their difference, and as before this is no bigger than the round-off error.

plot(xx, fnval(dpp,xx) - 3*subplus(xx-2).^2)
title('Error in Derivative of interpolant to ((x-2)_+)^3')

The second derivative of the truncated power is

f2′′(x)=6(x-2)+

A plot of the difference between this function and the second derivative of the interpolant to the original function shows that there are now jumps, but they are still within the round-off error.

ddpp = fnder(dpp);
plot(xx, fnval(ddpp,xx) - 6*subplus(xx-2))
title('Error in Second Derivative of Interpolant to ((x-2)_+)^3')

The integral of the truncated power is

F2(x)=((x-2)+)4/4.

A plot of the difference between this function and the integral of the interpolant to the original function again shows that the errors are within the round-off error.

ipp = fnint(pp);
plot(xx, fnval(ipp,xx) - subplus(xx-2).^4/4)
title('Error in Integral of Interpolant to ((x-2)_+)^3')

The CSAPE Command

Like csapi, the csape command provides a cubic spline interpolant to given data. However, it permits various additional end conditions. Its simplest version,

pp = csape(x,y)

uses the Lagrange end condition, which is a common alternative to the not-a-knot condition used by csapi. csape does not directly return values of the interpolant, but only its ppform.

For example, consider again interpolation to the function

f1(x)=((x-1)+)3,

which csapi fails to reproduce well. We plot the error of the not-a-knot interpolant returned by csapi (in black), along with the error of the interpolant obtained from csape (in red).

exact = subplus(xx-1).^3;
plot(xx, fnval(csapi(x,subplus(x-1).^3),xx) - exact,'k')
hold on
plot(xx, fnval(csape(x,subplus(x-1).^3),xx) - exact,'r')
title('Error in Not-a-Knot vs. Lagrange End Conditions')
legend({'Not-a-Knot' 'Lagrange'});
hold off

Figure contains an axes object. The axes object with title Error in Not-a-Knot vs. Lagrange End Conditions contains 2 objects of type line. These objects represent Not-a-Knot, Lagrange.

There is not much difference between the two interpolants in this case.

Other End Conditions: The ‘Natural’ Spline Interpolant

The csape command also provides ways to specify several other types of end conditions for an interpolating cubic spline. For example, the command

pp = csape(x,y,'variational')

uses the so-called ‘natural’ end conditions. This means that the second derivative is zero at the two extreme breaks.

This step shows how to apply ‘natural’ cubic spline interpolation to the function

f2(x)=((x-2)+)3,

pp = csape(x,subplus(x-2).^3,'second');
plot(xx, fnval(pp,xx) - subplus(xx-2).^3)
title('Error in ''Natural'' Spline Interpolation to ((x-2)_+)^3')

Note the large error near the right end. This is due to the fact that the ‘natural’ end conditions implicitly insist on having a zero second derivative there.

Other End Conditions: Prescribing Second Derivatives

We can also explicitly use the correct second derivatives to get a small error. First, we compute the correct second derivative values of the truncated power at the endpoints.

endcond = 6*subplus(x([1 end])-2);

pp = csape(x,[endcond(1) subplus(x-2).^3 endcond(2)], 'second');
plot(xx, fnval(pp,xx) - subplus(xx-2).^3,'r')
title(['Error in Spline Interpolation to ((x-1)_+)^3'; ...
       '  When Matching the 2nd Derivative at Ends  '])

Other End Conditions: Prescribing Slopes

csape also permits specification of endpoint slopes. This is the clamped (or, complete) cubic spline interpolant. The statement

pp = csape(x,[sl,y,sr],'clamped')

creates the cubic spline interpolant to the data (x, y) that also has slope sl at the leftmost data site and slope sr at the rightmost data site.

Other End Conditions: Mixed End Conditions

It is even possible to mix these conditions. For example, our much-exercised truncated power function

f1(x)=((x-1)+)3

has slope 0 at x=0 and second derivative 30 at x=6 (the last data site).

Therefore, by matching the slope at the left end and the curvature at the right, we expect no error in the resulting interpolant.

pp = csape(x, [0 subplus(x-1).^3 30], [1 2]);
plot(xx, fnval(pp,xx) - subplus(xx-1).^3)
title(['Error in Spline Interpolation to ((x-1)_+)^3'; ...
       '        with Mixed End Conditions.          '])

Other End Conditions: Periodic Conditions

x = (pi/2)*(-2:2);
y = [0 -1 0 1 0];
pp = csape(x,y, 'periodic' );
xx = linspace(-pi,pi,201);
plot(xx, sin(xx) - fnval(pp,xx), 'x')
title('Error in Periodic Cubic Spline Interpolation to sin(x)')

Figure contains an axes object. The axes object with title Error in Periodic Cubic Spline Interpolation to sin(x) contains an object of type line.

End Conditions Not Explicitly Covered by CSAPI or CSAPE

For example, suppose that you want to compute the cubic spline interpolant s to the data

x = 0:.25:3;
q = @(x) x.*(-1 + x.*(-1+x.*x/5));
y = q(x);

and enforce the condition

lambda(s) := a * (Ds)(e) + b * (D^2 s)(e) = c

on the first and second derivatives of s at the point e.

The data were generated from a quartic polynomial that happens to satisfy this side condition with specific parameters

e = x(1);
a = 2; b = -3; c = 4;

To construct the interpolant that satisfies this specific condition, we first construct the interpolant with the default end conditions

and the first derivative of its first polynomial piece.

dp1 = fnder(fnbrk(pp1,1));

In addition, we construct the cubic spline interpolant to zero data values, specifying that it have a slope of 1 at e,

pp0 = csape(x,[1,zeros(size(y)),0], [1,0]);

as well as constructing the first derivative of its first polynomial piece.

dp0 = fnder(fnbrk(pp0,1));

Then we compute lambda for both pp1 and pp0,

lam1 = a*fnval(dp1,e) + b*fnval(fnder(dp1),e);
lam0 = a*fnval(dp0,e) + b*fnval(fnder(dp0),e);

and construct the correct linear combination of pp1 and pp0 to get a cubic spline

s := pp1 + ((c - lambda(pp1))/lambda(pp0)) * pp0

that does satisfy the desired condition, as well as the default end condition at the right endpoint. We form this linear combination with the help of fncmb.

s = fncmb(pp0,(c-lam1)/lam0,pp1);

A plot of the interpolation error shows that s fits the quartic polynomial slightly better near e than the interpolant pp1 with the default conditions does.

xx = (-.3):.05:.7; yy = q(xx);
plot(xx, fnval(pp1,xx) - yy, 'x')
hold on
plot(xx, fnval(s,xx) - yy, 'o')
hold off
legend({'Default conditions' 'Nonstandard conditions'},'location','SE')

Figure contains an axes object. The axes object contains 2 objects of type line. These objects represent Default conditions, Nonstandard conditions.

If we want to enforce the condition

mu(s) := (D^3 s)(3) = 14.6

pp2 = csape(x,[0,zeros(size(y)),1],[0,1]);

Then we find the coefficients d0 and d2 in the linear combination

s := pp1 + d0*pp0 + d2*pp2

that solves the linear system

lambda(s) = c

mu(s) = 14.6

Note that both pp0 and pp2 vanish at all interpolation sites, hence s will match the given data for any choice of d0 and d2.

For amusement, we use the MATLAB® encoding facility to write a loop to compute lambda(pp_j) and mu(pp_j), for j=0:2.

dd = zeros(2,3);
for j=0:2
   J = num2str(j);
   eval(['dpp',J,'=fnder(pp',J,');']);
   eval(['ddpp',J,'=fnder(dpp',J,');']);
   eval(['dd(1,1+',J,')=a*fnval(dpp',J,',e)+b*fnval(ddpp',J,',e);']);
   eval(['dd(2,1+',J,')=fnval(fnder(ddpp',J,'),3);']);
end

Given the values of lambda and mu for pp0, pp1, and pp2, we then solve for the coefficients that define the correct linear combination.

d = dd(:,[1,3])([c;14.6]-dd(:,2));
s = fncmb(fncmb(pp0,d(1),pp2,d(2)),pp1);

xxx = 0:.05:3;
yyy = q(xxx);
plot(xxx, yyy - fnval(s,xxx),'x')
title('Error in Spline Interpolant to y = x*(-1 + x*(-1+x*x/5))')

Figure contains an axes object. The axes object with title Error in Spline Interpolant to y = x*(-1 + x*(-1+x*x/5)) contains an object of type line.

For reassurance, we compare this error with the one obtained in complete cubic spline interpolation to this function:

hold on
plot(xxx, yyy - fnval(csape(x,[-1,y,-7+(4/5)*27],'clamped'),xxx),'o')
hold off
legend({'Nonstandard conditions' 'Endslope conditions'})

The errors differ (and not by much) only near the end points, testifying to the fact that both pp0 and pp2 are sizable only near their respective end points.

As a final check, we verify that s satisfies the desired third derivative condition at 3.

Источник

#include <stdio.h>
#include <cstring>
#include <cmath>
#include <iostream>
 
class knot {
public:
    double x, f, f2;
    void Add(double arg, double func, double func2)
    {
        x = arg;
        f = func;
        f2 = func2;
    }
    knot() {}
};
 
class vector {
public:
    double* x;
    void Add(int m)
    {
        x = new double[m];
    }
    vector()
    {
    }
};
 
knot* KnotArray;
int n; // количество узлов интерполяции
double** Coef;
double* b;
 
//*************************************************************************/
//Решение системы уравнений с трехдиагональной матрицей
//*************************************************************************/
void SolveTriDiag(double** TDM, double* F)
{
    double* alph = new double[n - 1];
    double* beta = new double[n - 1];
 
    int i;
 
    alph[0] = -TDM[2][0] / TDM[1][0];
    beta[0] = F[0] / TDM[1][0];
 
    for (i = 1; i < n - 1; i++)
    {
        alph[i] = -TDM[2][i] / (TDM[1][i] + TDM[0][i] * alph[i - 1]);
        beta[i] = (F[i] - TDM[0][i] * beta[i - 1]) / (TDM[1][i] + TDM[0][i] * alph[i - 1]);
    }
    b[n - 1] = (F[n - 1] - TDM[0][n - 1] * beta[n - 2]) / (TDM[1][n - 1] + TDM[0][n - 1] * alph[n - 2]);
 
    for (i = n - 2; i > -1; i--)
    {
        b[i] = b[i + 1] * alph[i] + beta[i];
    }
}
 
//*************************************************************************/
//Построение таблицы коэффициентов кубического сплайна  y=f(x)
//*************************************************************************/
int BuildSpline()
{
    double* a = new double[n - 1];
    double* c = new double[n - 1];
    double* d = new double[n - 1];
    double* delta = new double[n - 1];
    double* h = new double[n - 1];
    double** TriDiagMatrix = new double* [3];
 
    b = new double[n];
 
    TriDiagMatrix[0] = new double[n];
    TriDiagMatrix[1] = new double[n];
    TriDiagMatrix[2] = new double[n];
 
    double* f = new double[n];
    double x3, xn;
    int i;
 
    if (n < 3)
        return -1;
 
    x3 = KnotArray[2].x - KnotArray[0].x;
    xn = KnotArray[n - 1].x - KnotArray[n - 3].x;
 
    for (i = 0; i < n - 1; i++)
    {
        a[i] = KnotArray[i].f;
        h[i] = KnotArray[i + 1].x - KnotArray[i].x;
        delta[i] = (KnotArray[i + 1].f - KnotArray[i].f) / h[i];
        TriDiagMatrix[0][i] = i > 0 ? h[i] : x3;
        f[i] = i > 0 ? 3 * (h[i] * delta[i - 1] + h[i - 1] * delta[i]) : 0;
    }
    TriDiagMatrix[1][0] = h[0];
    TriDiagMatrix[2][0] = h[0];
    for (int i = 1; i < n - 1; i++)
    {
        TriDiagMatrix[1][i] = 2 * (h[i] + h[i - 1]);
        TriDiagMatrix[2][i] = h[i];
    }
    TriDiagMatrix[1][n - 1] = h[n - 2];
    TriDiagMatrix[2][n - 1] = xn;
    TriDiagMatrix[0][n - 1] = h[n - 2];
 
 
    i = n - 1;
    f[0] = ((h[0] + 2 * x3) * h[1] * delta[0] + powf(h[0], 2) * delta[1]) / x3;
    f[n - 1] = (powf(h[i - 1], 2) * delta[i - 2] + (2 * xn + h[i - 1]) * h[i - 2] * delta[i - 1]) / xn;
 
    SolveTriDiag(TriDiagMatrix, f);
 
    /*Coef = new double* [4];
    Coef[0] = new double[n - 1];
    Coef[1] = new double[n - 1];
    Coef[2] = new double[n - 1];
    Coef[3] = new double[n - 1];
    */
    Coef = new double* [n - 1];
    for (int count = 0; count < n - 1; count++)
        Coef[count] = new double[4];
    int j;
 
    for (j = 0; j < n - 1; j++)
    {
        d[j] = (b[j + 1] + b[j] - 2 * delta[j]) / (h[j] * h[j]);
        c[j] = 2 * (delta[j] - b[j]) / h[j] - (b[j + 1] - delta[j]) / h[j];
 
        Coef[j][0] = a[j];
        Coef[j][1] = b[j];
        Coef[j][2] = c[j];
        Coef[j][3] = d[j];
    }
    /*for (j = 0; j < n - 1; j++)//вывод значений коэффициентов полиномов
    {
        for (int i = 0; i < 4; i++)
            printf("%lft", Coef[j][i]);
        printf("n");
    }*/
    return 1;
}
/*************************************************************************/
//Построение таблицы коэффициентов кубического сплайна z=f(x)
//*************************************************************************/
int BuildSpline_z()
{
    double* a = new double[n - 1];
    double* c = new double[n - 1];
    double* d = new double[n - 1];
    double* delta = new double[n - 1];
    double* h = new double[n - 1];
    double** TriDiagMatrix = new double* [3];
 
    b = new double[n];
 
    TriDiagMatrix[0] = new double[n];
    TriDiagMatrix[1] = new double[n];
    TriDiagMatrix[2] = new double[n];
 
    double* f = new double[n];
    double x3, xn;
    int i;
 
    if (n < 3)
        return -1;
 
    x3 = KnotArray[2].x - KnotArray[0].x;
    xn = KnotArray[n - 1].x - KnotArray[n - 3].x;
 
    for (i = 0; i < n - 1; i++)
    {
        a[i] = KnotArray[i].f2;
        h[i] = KnotArray[i + 1].x - KnotArray[i].x;
        delta[i] = (KnotArray[i + 1].f2 - KnotArray[i].f2) / h[i];
        TriDiagMatrix[0][i] = i > 0 ? h[i] : x3;
        f[i] = i > 0 ? 3 * (h[i] * delta[i - 1] + h[i - 1] * delta[i]) : 0;
    }
    TriDiagMatrix[1][0] = h[0];
    TriDiagMatrix[2][0] = h[0];
    for (int i = 1; i < n - 1; i++)
    {
        TriDiagMatrix[1][i] = 2 * (h[i] + h[i - 1]);
        TriDiagMatrix[2][i] = h[i];
    }
    TriDiagMatrix[1][n - 1] = h[n - 2];
    TriDiagMatrix[2][n - 1] = xn;
    TriDiagMatrix[0][n - 1] = h[n - 2];
 
 
    i = n - 1;
    f[0] = ((h[0] + 2 * x3) * h[1] * delta[0] + powf(h[0], 2) * delta[1]) / x3;
    f[n - 1] = (powf(h[i - 1], 2) * delta[i - 2] + (2 * xn + h[i - 1]) * h[i - 2] * delta[i - 1]) / xn;
 
    SolveTriDiag(TriDiagMatrix, f);
 
    /*Coef = new double* [4];
    Coef[0] = new double[n - 1];
    Coef[1] = new double[n - 1];
    Coef[2] = new double[n - 1];
    Coef[3] = new double[n - 1];*/
 
    Coef = new double* [n - 1];
    for (int count = 0; count < n - 1; count++)
        Coef[count] = new double[4];
 
 
    int j;
    for (j = 0; j < n - 1; j++)
    {
        d[j] = (b[j + 1] + b[j] - 2 * delta[j]) / (h[j] * h[j]);
        c[j] = 2 * (delta[j] - b[j]) / h[j] - (b[j + 1] - delta[j]) / h[j];
 
        Coef[j][0] = a[j];
        Coef[j][1] = b[j];
        Coef[j][2] = c[j];
        Coef[j][3] = d[j];
    }
 
}
//*************************************************************************/
//Подсчет значения интерполянты в заданной точке
//*************************************************************************/
double Interpolate(double x)
{
    //double result;
    int i = 0;
 
    while (KnotArray[i].x < x)
        i++;
    i--;
    return Coef[i][0] + Coef[i][1] * (x - KnotArray[i].x) + Coef[i][2] * powf((x - KnotArray[i].x), 2) + Coef[i][3] * powf((x - KnotArray[i].x), 3);
 
}
 
//*************************************************************************/
//Загрузка данных
//*************************************************************************/
int Load_Data()
{
    printf("Input filename with datan");
    char FileName[20];
    int i = 0;
    FILE* File;
    scanf_s("%s", &FileName, 20);
    if (!fopen_s(&File, FileName, "r"))
    {
        printf("file is openedn");
    }
    else
    {
        printf("%s: file doesn't existn", FileName);
        return -1;
    }
    double x, f, f2;
    fscanf_s(File, "%d", &n);
    KnotArray = new knot[n + 2];
    while (!feof(File))
    {
        fscanf_s(File, "%lf%lf%lf", &x, &f, &f2);
        KnotArray[i].Add(x, f, f2);
        i++;
        if (i == n)
            return 1;
    }
    fclose(File);
    return 1;
}
 
int main()
{
    FILE* file_out;
 
    int N;//количество значений функции-интерполянты
    //double x = 0;
    printf("Input the number of interpolant values:n");
    scanf_s("%d", &N);
    double* init = new double[N + 1];
    double* fun = new double[N + 1];
    double* fun_z = new double[N + 1];
 
    if (Load_Data() != -1)
    {
        for (int i = 0; i <= N; i++)
        {
            init[i] = KnotArray[0].x + (KnotArray[n - 1].x - KnotArray[0].x) / N * i;
 
        }
        BuildSpline();
        fun[0] = KnotArray[0].f;
        fun_z[0] = KnotArray[0].f2;
        fun[N] = KnotArray[n - 1].f;
        fun_z[N] = KnotArray[n - 1].f2;
        for (int i = 1; i < N; i++)
        {
            fun[i] = Interpolate(init[i]);
            printf("i=%d f=%lfn", i, fun[i]);
        }
        BuildSpline_z();
        for (int i = 1; i < N; i++)
        {
            fun_z[i] = Interpolate(init[i]);
        }
        if (fopen_s(&file_out, "D:interp.dat", "wb"))
            printf("File could not be openedn");
        else
        {
            for (int i = 0; i <= N; i++)
            {
                printf("x:  %lfty_interp:   %lftz_interp:   %lfn", init[i], fun[i], fun_z[i]);
                fwrite(&init[i], sizeof(double), 1, file_out);
                fwrite(&fun[i], sizeof(double), 1, file_out);
                fwrite(&fun_z[i], sizeof(double), 1, file_out);
            }
        }
    }
    delete[] init;
    delete[] fun;
    delete[] fun_z;
    for (int count = 0; count < n - 1; count++)
        delete[] Coef[count];
    system("pause");
    return 0;
}

Источник

Правила форума

В этом разделе нельзя создавать новые темы.

Если Вы хотите задать новый вопрос, то не дописывайте

его в существующую тему, а создайте новую в корневом разделе «Помогите решить/разобраться (М)».

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

Не ищите на этом форуме халяву

, правила запрещают участникам публиковать готовые решения стандартных учебных задач. Автор вопроса обязан привести свои попытки решения

и указать конкретные затруднения.

Обязательно просмотрите тему

Правила данного раздела, иначе Ваша тема может быть удалена

или перемещена в Карантин, а Вы так и не узнаете, почему.

	Погрешность интерполяции сплайнами и оптимальная сетка 18.10.2007, 15:57
04/01/07 90	Всем привет! Вопрос к тем, кто изучал погрешность сплайн-интерполяции. Известно, что оценить ее сверху можно с использованием производных высших степеней интерполируемой функции. Используя эти сведения я попытался состряпать алгоритм формирования оптимальной сетки измерений. Но вскоре понял, что такой подход — тупик, поскольку производные (тем более высоких степеней) известны мне с погрешностью ~ +- 100% . Решил попробовать моделировать процедуру интерполяции для типовой (средне-обычной) функции и подбирать сетку таким образом, чтоб погрешность не превысила допустимой. Тут своии проблемы, связанные с взаимным влиянием интервалов. Т.е. при изменении шага измерения на одном из интервалов, меняется погрешность на другом и наоборот. В принципе, очень грубо, тут можно прикинуть в каких точках мерять, но слишком халтурно и «на глаз». Если кто-нибудь занимался подобной деятельностью, посоветуйте, пожалуйста, более-менее обоснованный способ расчета оптимальной сетки измерений. Заранее благодарен.

worm2

Re: Погрешность интерполяции сплайнами и оптимальная сетка

18.10.2007, 16:29

Заслуженный участник

01/08/06
2928
Уфа

oliva писал(а):

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

Что значит «высших степеней»? 3-я это высшая или как?

Для сплайн-интерполяции

погрешность зависит от порядка сплайна, но не от количества точек. Например, для кубического сплайна на равномерной сетке:

$|R_n(x)| leqslant frac{5}{2} h^3 maxlimits_{y in [a, b]} |f^{(3)}(y)|$ , где h — шаг сетки.

Вас такая погрешность не устраивает?

Может быть, вы имеете в виду полиномиальную интерполяцию (единым многочленом на всём отрезке) ?.

Добавлено спустя 2 минуты 49 секунд:

Или у Вас нет вообще никаких оценок для максимума модуля третьей производной?

Профиль

oliva	Re: Погрешность интерполяции сплайнами и оптимальная сетка 18.10.2007, 16:46
04/01/07 90	worm2 писал(а): для кубического сплайна на равномерной сетке: $\|R_n(x)\| leqslant frac{5}{2} h^3 maxlimits_{y in [a, b]} \|f^{(3)}(y)\|$ , где h — шаг сетки. Вообще-то у меня другая формула для погрешности кубического сплайна. Но спорить не буду. Я так понял, что эта погрешность определяется свойствами (прерывностью, дифференциируемостью и т.п.) функции f. worm2 писал(а): Или у Вас нет вообще никаких оценок для максимума модуля третьей производной? Оценить максимум модуля я могу, но приблизительно. При этом получаю низкую точность определения шага измерения. Особенно это заметно в местах резкого изменения свойств f Добавлено спустя 3 минуты 28 секунд: worm2 писал(а): погрешность зависит от порядка сплайна, но не от количества точек. Например, для кубического сплайна на равномерной сетке: $\|R_n(x)\| leqslant frac{5}{2} h^3 maxlimits_{y in [a, b]} \|f^{(3)}(y)\|$ , где h — шаг Шаг и количество точек я рассматриваю как взаимоопределяющие факторы.

olga_helga	18.10.2007, 17:28
28/09/07 86	Стоп, я так понимаю что по сути это все дело сводится к нахождению новой таблицы значений, по которой будет построен более точный сплайн.Так? Добавлено спустя 4 минуты 17 секунд: А если так,то по сути это есть «задача зглаживания эксперементальных данных», основная идея которой построение новой таблицы данных, путем получения нового сглаженного значения сведением какой либо зависимости м/у и к линейной

oliva	18.10.2007, 17:29
04/01/07 90	olga_helga писал(а): Стоп, я так понимаю что по сути это все дело сводится к нахождению новой таблицы значений, по которой будет построен более точный сплайн.Так? Речь идет об интерполяции данных, полученных в процессе градуировки датчиков (в данном случае температуры). Характеристика этих датчиков в рабочем диапазоне меняется от почти линейной до сильно нелинейной. Нужно составить оптимальную сетку измерений, которая гарантирует нужную точность интерполяции при минимальном количестве измерений. Для этого я взял типовую характеристику и пытаюсь ее исследовать.

olga_helga	18.10.2007, 17:44
28/09/07 86	Вооо!!!! Тык вот, для зглаживания экспериментальных данных используют апроксимацию. Используя старую таблицу данных $[ left( {x_i ,y_i } right) ]$ , получают новую $[ left( {x_i ,overline {y_i } } right) ]$ , где $[ {overline {y_i } } ]$ -сглаженное значение. Так вот коим образом получают новую таблицу. 1)если зависимость линейная, то $[ overline {y_i } = frac{{y_{i - 1} + y_i + y_{i + 1} }} {3} ]$ 2)если зависимость степенная $[ y = ae^{bx} ]$ ,,то вводя замену переменных , получаем ,а $[ left( {ln x_i ,ln y_i } right) ]$ -новый узел таблицы 3)если зависимость $[ y = frac{a} {{b + x}} ]$ , то замена , -линейная зависимость, а $[ left( {x_i cdot y_i ,y_i } right) ]$ -новая таблица. Ну вот что-то вроде етого.

oliva	18.10.2007, 17:52
04/01/07 90	olga_helga писал(а): для зглаживания экспериментальных данных используют апроксимацию. Да не нужно данные сглаживать Аппроксимацию со сглаживанием применяют когда данные получены с погрешностью и им не доверяют. Я решаю задачу интерполяции и измеренным точкам доверяю больше чем каким-либо расчетным. Интерполяция — нахождение кривой, которая ПРОХОДИТ через экспериментальные точки.

olga_helga	18.10.2007, 18:14
28/09/07 86

worm2

18.10.2007, 18:17

Заслуженный участник

01/08/06
2928
Уфа

Правильно ли я понял задачу:

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

Интересная задача. Можно попробовать решать тупым перебором на компьютере с каким-то шагом.

Добавлено спустя 2 минуты 37 секунд:

Я понимаю, что Ваша задача намного сложнее, но с этой упрощённой формулировки можно было бы начать.

Профиль

olga_helga	18.10.2007, 18:21
28/09/07 86	Многочлен имеет n действительных корней,выражаемых формулой $[ x_i = cos frac{{(2i + 1)pi }} {{2n}},i = 0,1,...,n - 1 ]$ Добавлено спустя 3 минуты 14 секунд: В роде есть еще формулка для $[ x_i^{} in [a,b] ]$ ,что то типа что узлы интерполяции выбираются как $[ x_i = frac{1} {2}left( {(b - a)cos frac{{(2i + 1)pi }} {{2n + 2}} + b + a} right),i = 0,1,...,n ]$ .Воть.Помогло?

Brukvalub

18.10.2007, 18:22

Заслуженный участник

01/03/06
13626
Москва

olga_helga писал(а):

Помоему есть такая оценка погрешности интерполяции $[ left| {f(x) - L_n (x)} right| leqslant frac{{mathop {max }limits_{x in [a,b]} }} {{(n + 1)!}}left| {prodlimits_{i = 0}^n {(x - xi)} } right| ]$

Это неверная оценка. Необходимо еще потребовать существование $[f^{(n + 1)} (x)]$ и умножить правую часть на $[left| {f^{(n + 1)} (xi )} right|]$

Профиль

olga_helga	18.10.2007, 18:27
28/09/07 86	Извиняюсь!Я ошибочку сделала! В оригинале оценка для интерполирования n+1 раз непрерывно дифференцируемой ф-ции на отрезке [a,b] ф-ции f(x) многочленом Логранжа такая: $[ left\| {f(x) - L_n (x)} right\| leqslant frac{{mathop {max left\| {f^{(n + 1)} (x)} right\|}limits_{x in [a,b]} }} {{(n + 1)!}}left\| {prodlimits_{i = 0}^n {(x - x_i )} } right\| ]$ .Вроде так.

Brukvalub

18.10.2007, 18:29

Заслуженный участник

01/03/06
13626
Москва

Профиль

olga_helga	to worm2 18.10.2007, 18:31
28/09/07 86	Попробуйте последней из написанных формул воспользоваться.

oliva	Re: to worm2 18.10.2007, 19:36
04/01/07 90	olga_helga писал(а): Попробуйте последней из написанных формул воспользоваться. Я ж начал с того, что не могу с достаточной точностью знать производную выше первой степени, а в Вашей формуле присутствует производная степени (n+1). Что-же касается полиномов Чебышева — не все понял, но интересно. Сообщите пожалуйста источник, чтоб я вник. Добавлено спустя 5 минут 24 секунды: worm2 писал(а): Можно попробовать решать тупым перебором на компьютере с каким-то шагом. Над этим и «пыхтю» Проблема в том, что интервал измерения переменный и влияет на соседние. Т.е. при попытке подобрать один интервал (если подходить строго) нужно пересчитывать остальные, что усложняет, запутывает и ставит под сомнение результативность алгоритма.

Модераторы: Модераторы Математики, Супермодераторы

Источник

Ошибка интерполяции кубическими сплайнами

Материал из MachineLearning.

Содержание

Введение

Постановка математической задачи

Изложение метода

Метод прогонки

Пример: интерполирование неизвестной функции

Ошибка интерполяции

Пример: интерполяция синуса

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

См. также

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

3. Локальный сплайн.

4. Глобальные способы построения кубических сплайнов.

5. Погрешность приближения кубическими сплайнами.

The CSAPI Command

How to Check the Output of CSAPI

Example: The Truncated Power Function

The Interpolation Error

A Truncated Power That Cannot be Reproduced

Using the ppform Instead of Values

Example: Differentiating and Integrating the Interpolant

The CSAPE Command

Other End Conditions: The ‘Natural’ Spline Interpolant

Other End Conditions: Prescribing Second Derivatives

Other End Conditions: Prescribing Slopes

Other End Conditions: Mixed End Conditions

Other End Conditions: Periodic Conditions

End Conditions Not Explicitly Covered by CSAPI or CSAPE

The CSAPI Command

How to Check the Output of CSAPI

Example: The Truncated Power Function

The Interpolation Error

A Truncated Power That Cannot be Reproduced

Using the ppform Instead of Values

Example: Differentiating and Integrating the Interpolant

The CSAPE Command

Other End Conditions: The ‘Natural’ Spline Interpolant

Other End Conditions: Prescribing Second Derivatives

Other End Conditions: Prescribing Slopes

Other End Conditions: Mixed End Conditions

Other End Conditions: Periodic Conditions

End Conditions Not Explicitly Covered by CSAPI or CSAPE

Правила форума

А вот еще интересные материалы: