1.      »нтерпол€ци€ сплайнами

 

ѕри интерполировании многочленами с увеличением числа узлов интерпол€ции увеличиваетс€ и степень интерпол€ционного многочлена. ѕри увеличении узлов интерпол€ции точность вычислений растет, но при этом величины производных увеличиваютс€ еще быстрее, и в точках отличных от узлов интерпол€ции возрастает ошибка вычислени€. ¬ последние годы по€вилась иде€ интерполировать сплайнами. —плайн в переводе на русский €зык рейка. ћожно интерполировать отрезками пр€мых, соедин€ющих узлы, но при этом интерпол€нтом €вл€етс€ недостаточно гладка€ крива€ (имеет изломы). ѕо€вилась иде€ интерполировать кусками полиномов невысокой степени (второй или третьей). “ак чтобы получилась непрерывна€ гладка€ крива€, производные на стыках кусков должны совпадать.

Ѕудем говорить о приближении функции f(x) на отрезке [a, b]. ѕусть на отрезке заданы узлы интерпол€ции a= x0< x1 <Е < xn = b, в которых известны значени€ функции f(x) . «адача интерполировани€ сплайнами формулируетс€ следующим образом:

Ќайти функцию Sm(x), удовлетвор€ющую требовани€м:

1) Sm(x) на каждом из отрезков [xk-1, xk ], k = 1,Е, n, €вл€етс€ многочленом Pmk(x) m - й степени:

Pmk(x) = amk xm + am-1k xm-1+ Е+ a1k x + a0k ;

2) в узлах интерпол€ции xk имеет место равенство

Sm(xk) = f(xk); k = 0, 1, Е, n;

3) Sm(x) на [a, b] имеет непрерывные производные до m - 1 -го пор€дка включительно, т.е. выполн€ютс€ равенства

Pmk(xk).

‘ункци€ называетс€ сплайном m - го пор€дка.

¬ нашем распор€жении (m +1)n неизвестных коэффициентов ami искомых многочленов Pmk(x). ”слови€ 2) и 3) дают нам 2 n + (m - 1)(n -1) = (m - 1)n + n +1 уравнений дл€ определени€ этих коэффициентов. Ќедостающие m - 1 уравнений дописывают, накладыва€ краевые услови€. ≈сли m >1 , то дл€ единственности Sm(x) следует задавать ещеусловий, которые задаютс€ на концах отрезка либо произвольно, либо из дополнительной информации о поведении функции .

ѕри m =1 получаем интерпол€цию отрезками (метод ломаных). ѕри этом Sm(x) равномерно сходитс€ к непрерывной функции f(xk) на [a, b], если при n Ѓ ¥.

ѕри m =2получаем интерпол€цию дугами парабол (квадратичный сплайн).

Ќаиболее часто используетс€ кусочно-кубическа€ интерпол€ци€, т.е. интерпол€ци€ кубическими сплайнами.

«аметим, что при увеличении пор€дка m сплайна скорость сходимости повышаетс€ с увеличением пор€дка сплайна и гладкости f(x).

–ассмотрим построение кубического сплайна S3(x) на отрезке [a, b], на котором заданы значени€ функции f(xk) в узлах x0, x1,Е , xn. ѕо определению S3(x) на отрезке [xk-1, xk ], k = 1,Е, n, €вл€етс€ кубическим многочленом. “огда имеем

1) P3k(x) = a3k x3 + a2k x2+ a1k x + a0k ; x Î [xk-1, xk ]; k = 1,Е, n;

2) S3(xk) = f(xk); k = 0, 1, Е, n; т.е. должны выполн€тс€ услови€:

P3k(xk) = P3k+1(xk) = f(xk), k = 1, Е, n-1;

P31(x0) = f(x0), P3n(xn) = f(xn);

(всего 2n соотношений), или учитыва€ вид P3k(x):

a3k xk3 + a2k xk2+ a1k xk + a0k= f(xk), k = 1, Е, n-1;

a3k+1 xk3 + a2k+1 xk2+ a1k +1xk + a0k+1= f(xk), k = 1, Е, n-1;

a31 xk3 + a21 xk2+ a11 xk + a01= f(x0),

a3n xk3 + a2n xk2+ a1n xk + a0n= f(xn);

3)

т.е

3a3k xk2 + 2a2k xk+ a1k= 3a3k+1 xk2 + 2a2k+1 xk+ a1k +1, k = 1, Е, n-1;

6a3k xk + 2a2k ††= 6a3k+1 xk + 2a2k+1 , k = 1, Е, n-1.

Ќаложим еще два краевые услови€:

или

6a31x0 + 2a21 ††= 0,

6a3n x n + 2a2n = 0.

ѕолучим —Ћј”, дл€ решени€ которой применим следующий прием. ќбозначим через Mk вторую производную S3(x) во внутренней точке xk, т.е.

Mk = S3"( x k); k = 1, Е, n-1.

”читыва€ непрерывность вторых производных Sm(x) в точках x k ; k = 1, Е, n-1, дл€ Mk имеем

Mk =6a3k x k + 2a2k ,

и

Mk =6a3k +1 x k + 2a2k+1.

“ак как S3"( x ) на каждом из интервалов (x k-1, x k) €вл€етс€ линейной функцией от x:

S3"( x ) = 6a3k x + 2a2k , xÎ(x k-1, x k),

и, следовательно, S3"( x ) можетбыть представлена через значени€ Mk -1 и Mkна концах

или

.

“огда

.

ƒважды интегриру€ полученное выражение дл€ S3"( x ) и учитыва€, что S3(x k-1) = f(xk-1), S3(x k) = f(xk), дл€ S3 (x) при x Î[x k-1, x k] получим следующее выражение:

†††††††††††††††††††††††††††††††††††††††††††††† (5.1)

ƒифференциру€ S3(x) и учитыва€ непрерывность первых производных в точках x k ; k = 1, Е, n-1, т.е. услови€ получим систему из n-1 уравнений относительно M0, M1, Е, Mn:

”читыва€, что M0 = 0, Mn = 0, приходим к системеиз n-1 уравнений с n-1 неизвестными M1, M2 Е, Mn-1, матрица которой симметрична и трехдиагональна.

–еша€, полученную —Ћј”, методом √аусса или методом прогонки, получим значени€ M1, M2 Е, Mn-1. ѕодставив эти значени€ в (5.1) определим многочлены P3k(x), k = 1, Е, n, т.е. определим сплайн S3(x).

“еорема 4.1. »з всех функций, удовлетвор€ющих услови€м 2)-3) сплайн о пор€дка дает минимум дл€ функционала

.

јлгоритм интерпол€ции функции кубическими сплайнами.

¬вод: ”злы интерпол€ции X[i], Y[i]' i = 0,1,Е, n. “очка x.

¬ывод: ¬ычислить c:= S3(x).

÷икл по j:= 1Е n выполнить;

†††††† еслиi =1, то A[i]:= X[i]/6,B[i]:= (X[i+1])/3;C[i]:= (X[i+1]-X[i])/6;D[i]:= (Y[i+1]-Y[i])/ (X[i+1]-X[i]) - Y[i]/ Y[i];

иначе еслиi =n, то A[i] := (X[i]-X[i-1])/6,B[i]:= (-X[i-1])/3;C[i]:=(-X[i])/6; D[i]:= Y[i]/ X[i] - (Y[i]-Y[i-1])/ (X[i]-X[i-1]);

иначе A[i] := (X[i]-X[i-1])/6; B[i]:= (X[i+1] -X[i-1])/3; C[i]:= (X[i+1]-X[i])/6; D[i]:= (Y[i+1]-Y[i])/ (X[i+1]-X[i]) -- (Y[i]-Y[i-1])/ (X[i]-X[i-1]);

конец цикла по i;

÷икл по j:= 1Е n выполнить;

если j=1 то P[j]:= -C[j]/B[j]; L[j]:= -D[j]/B[j];

иначе P[j]:= -C[j]/(B[j]+A[j]*P[j-1]; L[j]:= (D[j]-A[j]*L[j-1])/(B[j]+A[j]*P[j-1]);

конец цикла по j;

M[0]:=0; M[n]:=0;

÷икл по j:= n-1Е 1 выполнить;

если j=n-1, то M[j]:= L[j];

иначе M[j]:= P[j]*M[j+1]+l[j];

конец цикла по j;

÷икл по k:= 1Е n выполнить;

если x>=X[k-1] и x<=X[k] , то

††††††††† если k=1, то с:= M[k-1]*( X[k]-x)* ( X[k]-x)* ( X[k]-x)/ ((X[k]-X[k-1])*6);

c:=c+ M[k]*( x- X[k-1])* (x- X[k-1])* (x- X[k-1])/ ((X[k]-X[k-1])*6);

c:=c+ (Y[k-1]-M[k-1]*( X[k]- X[k-1])* (x- X[k])* (X[k]-x))/ ((X[k]-X[k-1])*6);

c:=c+ (Y[k]-M[k]*( X[k]- X[k-1])* (x- X[k-1])* (x-X[k-1]))/ ((X[k]-X[k-1])*6); конец

конец цикла по k;

 

ѕример. ѕо таблице значений функции вычислить f(10)

x

2

5

7

8

12

f(x)

12

33

-12

2

20

 

 

xi

f(xi)

ai

bi

ci

di

pi

li

Mi

xi

2

5

7

8

12

12

33

-12

2

3

3

2

1

4

7

-22,5

14

0,25

5

3

5

-5,9

12,2

9,3

6

7

3,0

0,4

10

-0,3

f(10) = 120.