Введение: зачем и как нужно интерполировать

Первая проблема, с которой сталкиваются математики, работающие в приложениях, состоит в том, что предлагамые им исходные данные устроены совсем не так, как их учили на мехмате. Нет ни функций, ни формул, есть просто результаты тех или иных измерений, представляющие собой числа, организованные с той или иной степенью разумности. В частности, при решении прикладных геометрических задач, мы получаем на входе координаты точек в трехмерном пространстве, которые нам самим предстоит организовать в кривые или поверхности. Поэтому каждый раз приходится начинать с решения задачи ИНТЕРПОЛЯЦИИ, которая как раз и состоит в восстановлении кривой или поверхности по данному набору точек. При этом, полученная в результате кривая или поверхность должны обладать разумными математическими свойствами (например, быть дважды непрерывно дифференцируемой), не отходить далеко от исходных точек, а также выглядеть достаточно эстетично. Ясно, что от того, насколько хорошо мы справимся с этой задачей, зависит точность и эффективность всех наших дальнейших рассчетов. Тем самым, интерполяция имеет первостепенное значение в приложениях.

Рассмотрим простой пример. На первом курсе нам рассказывали как интерполировать график функции одного вещественного переменного с помощью интерполяционного многочлена. Для этого, напомним, надо выбрать на графике (n+1)  точку и построить единственный многочлен степени n, график которого через эти точки проходит.  При этом, качество полученной интерполяции существенно зависит от расположения выбранных точек. Если выбрать точки на графике, например, равномерно, то отклонение итогового  многочлена от приближаемой функции может расти с ростом степени полинома.

DynamicModule[{RungePlot,r,Runge,x},
Manipulate[
RungePlot[n],
{{n,5,"количество точек"},5,50,5,Appearance->"Labeled"},SaveDefinitions->True],Initialization:> (
r[x_]:=(1/(1 + 25 (x^2)));
Runge[n_]:=Table[ {x, r[x]}, {x, -1., 1., 2./(n-1)}];
RungePlot[n_]:=Module[{data,interp},
interp =Expand[InterpolatingPolynomial[Runge[n],x]];
Column[{Plot[{1/(1+25x^2),interp},{x,-1,1},ImageSize->{450},Epilog->{PointSize[Medium],Red,Point/@Runge[n]},PlotStyle->{Orange,Blue},PlotRange->{-1,1},Ticks->{{-1,1},{-1,1}},ImageSize->400,PlotRangePadding->.2],
ListLinePlot[Table[{x,Abs[N[N[interp] - r[x]]]},{x,-1,1,.01}],ImageSize->400,PlotRange->All,Ticks->{{-1,1},Automatic},PlotLabel->"Ошибка интеполяции",PlotStyle->Black,AspectRatio->.25,ImageSize->450]},Center]
])]

Проиллюстрованный здесь эффект носит название феномена Рунге. Преодолеть эту трудность можно с помощью более хитрого выбора точек на графике. Так например, если выбирать точки известным методом Чебышева, то для аналитических функций погрешность убывает с ростом степени многочлена  как геометрическая прогрессия. Другой широко применяемый метод состоит в использовании кусочно-полиномиальных интерполяционных функций – так называемых сплайнов.

В дословном переводе слово "spline" означает планку или рейку, используемую в строительных или морских конструкциях. Впервые понятие сплайна как кусочно-полиномиальной функции было введено  Исаакаом Шёнбергом (1903-1990, Румыния, США) в 1943-1946 годах при решении прикладных задач для американских ВВС, а затем  – и в теоретических исследованиях. Широкую известность и применение сплайны нашли после того, как в 60-е годы прошлого века французские математики-прикладники Пьер Безье, работавший на Рено,  и Поль де Кастелье, работавший на Ситроен, независимо применили их для создания первых САПР в автомобилестроении. Конечно, мы не знаем, как выглядела их программа, ну суть ее вероятно сводилась к следующей демо-версии.

DynamicModule[{color,carbody,degrees,mesh},
Manipulate[
Graphics3D[{
color,Specularity[White,50],
BSplineSurface[carbody,SplineDegree->Floor[degrees]],
If[mesh,
{Dashed,Blue,Line[carbody],Line[Transpose[carbody]],Red,PointSize[Medium],Point/@carbody},{}]},
Boxed->False,Lighting->"Neutral",RotationAction->"Clip",ViewPoint->{-2,-2,2}],
Row[{Control[{{degrees,{1,1}},{1,1},{14,7},ImageSize->{140,70}}],Dynamic[Floor[degrees]]},Spacer[10]],
Row[{Control[{color,ColorData["HTML","GoldenRod"]}],Control[{mesh,{False,True}}]},Spacer[10]]
],Initialization:>(carbody={
{{-.2,-2.25,0},{-.2,-2.25,0},{-.2,-2,0},{-.2,0,0},{-.2,2,0},{-.2,2.25,0},{-.2,2.25,0}},
{{0,-2.25,0},{0,-2.25,1.75},{0,-2,1.75},{0,0,1.75},{0,2,1.75},{0,2.25,1.75},{0,2.25,0}},
{{1,-2.25,0},{1,-2.25,1.75},{1,-2,1.75},{1,0,1.75},{1,2,1.75},{1,2.25,1.75},{1,2.25,0}},
{{1.01,-2.25,1.5},{1.01,-2.25,1.75},{1.01,-2,1.75},{1.01,0,1.75},{1.01,2,1.75},{1.01,2.25,1.75},{1.01,2.25,1.5}},
{{2.99,-2.25,1.5},{2.99,-2.25,1.75},{2.99,-2,1.75},{2.99,0,1.75},{2.99,2,1.75},{2.99,2.25,1.75},{2.99,2.25,1.5}},
{{3,-2.25,0},{3,-2.25,1.75},{3,-2,1.75},{3,0,1.75},{3,2,1.75},{3,2.25,1.75},{3,2.25,0}},
{{4,-2.25,0},{4,-2.25,1.75},{4,-2,1.75},{4,0,1.75},{4,2,1.75},{4,2.25,1.75},{4,2.25,0}},
{{5,-2.25,0},{5,-2.25,1.75},{5,-2,3.5},{5,0,3.5},{5,2,3.5},{5,2.25,1.75},{5,2.25,0}},
{{8,-2.25,0},{8,-2.25,1.75},{8,-2,3.5},{8,0,3.5},{8,2,3.5},{8,2.25,1.75},{8,2.25,0}},
{{8.01,-2.25,1.5},{8.01,-2.25,1.75},{8.01,-2,3.5},{8.01,0,3.5},{8.01,2,3.5},{8.01,2.25,1.75},{8.01,2.25,1.5}},
{{9,-2.25,1.5},{9,-2.25,2},{9,-2,2},{9,0,2},{9,2,2},{9,2.25,2},{9,2.25,1.5}},
{{9.99,-2.25,1.5},{9.99,-2.25,2},{9.99,-2,2},{9.99,0,2},{9.99,2,2},{9.99,2.25,2},{9.99,2.25,1.5}},
{{10,-2.25,0},{10,-2.25,2},{10,-2,2},{10,0,2},{10,2,2},{10,2.25,2},{10,2.25,0}},
{{11,-2.25,0},{11,-2.25,2},{11,-2,2},{11,0,2},{11,2,2},{11,2.25,2},{11,2.25,0}},
{{11.01,-2.25,0},{11.01,-2.25,0},{11.01,-2,0},{11.01,0,0},{11.01,2,0},{11.01,2.25,0},{11.01,2.25,0}}
})
]

Инженер-проектировщик задает контуры будущего автомобиля набором точек, которые отражают лишь основные функциональные особености (наличие четырёх ниш колеса, кабины, капота и пр.) Затем программа строит по эти точкам поверхность-сплайн, свойствами которого (в нашем примере – степенями многочленов-кусков по координатам) можно управлять. После того как  сплайн выбран, начинается работа по деталировке.

На первых нескольких занятиях этого семестра мы познакомимся с базовыми понятиями современной теории сплайнов. Разумеется, у нас нет времени сколько-нибудь углублятся в подробности. Современная теория сплайнов представляет собой хорошо разработанную, разветвленную и сложную науку. Заинтересовавшиеся слушатели приглашаются на спецкурс кафедры Дифференциальной геометрии и приложений "Компьютерная геометрия", который читают Д.П.Ильютко и Г.В.Носовский. См. также книгу Н.Н.Голованов, Д.П.Ильютко, Г.В.Носовский, А.Т.Фоменко, "Компьютерная геометрия''.

Мы начнем с одномерного случая.

.1. СПЛАЙНЫ

Определение. СПЛАЙНОМ называется кривая, построенная по совокупности точек пространства (Subscript[p, 0]),… ,(Subscript[p, n]), которые называются ОПОРНЫМИ ТОЧКАМИ сплайна. Требуется, чтобы сплайн r=r(t) В ТОЧНОСТИ проходил через свои опорные точки (Subscript[p, 0]),… ,(Subscript[p, n]) при заданных значениях параметра (Subscript[t, 0])≤ … ≤ (Subscript[t, n]) :  r((Subscript[t, i]))=(Subscript[p, i]), 0 ≤ i≤ n. Эти значения параметра (Subscript[t, i]) (0 ≤ i≤ n) называются УЗЛАМИ сплайна.

ПРОГРАММА 1. Ломаная, проходящая через заданные опорные точки  (Subscript[p, 0]),... ,(Subscript[p, n]).

DynamicModule[{pts},Manipulate[Show[Graphics[{Blue,Line[{pts}]}]],{{pts,{{1,1},{2,2.3},{3,2.7},{4,3},{5,4}}},Locator}]]

Рассмотрим основные примеры сплайнов.

.1.1. СПЛАЙН ЛАГРАНЖА

СПЛАЙН ЛАГРАНЖА аналогичен обычному интерполяционному многочлену Лагранжа. Это полиномиальная по параметру (и, следовательно бесконечно гладкая на всей области определения [(Subscript[t, 0]),(Subscript[t, n])]) кривая, которая задается формулой

r(t)=(∑_(i=0)^n)(Subscript[L, i])(t)(Subscript[p, i ]), t∈[(Subscript[t, 0]),(Subscript[t, n])],

где коэффициенты (Subscript[L, i])(t) зависят от узлов сплайна и задаются формулой:

(Subscript[L, i])(t)=(((t-(Subscript[t, 0]))…(t-(Subscript[t, i-1]))(t-(Subscript[t, i+1]))…(t-(Subscript[t, n])))/(((Subscript[t, i])-(Subscript[t, 0]))…((Subscript[t, i])-(Subscript[t, i-1]))((Subscript[t, i])-(Subscript[t, i+1]))…((Subscript[t, i])-(Subscript[t, n]))))=(∏_(j=0, j ≠ i)^n)(t-(Subscript[t, j]))/(∏_(j=0, j ≠ i)^n)((Subscript[t, i])-(Subscript[t, j])).

Очевидно, что для (Subscript[L, i])(t) выполняются равенства (Subscript[L, i])((Subscript[t, j]))=(Subscript[δ, ij]) , где (Subscript[δ, ij])  –  символ Кронекера.

Главный недостаток сплайнов Лагранжа состоят в том, что между двумя соседними опорными точками они могут образовывать петли. Такие петли, не обусловленные расположенем опорных точек, называются НЕОБОСНОВАННЫМИ ПЕТЛЯМИ и являются существенным недостатком сплайна, не позволяющим применять данный сплайн для большинства приложений компьютерной геометрии.

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

Эти недостатки – цена за бесконечную гладкость сплайна НА ВСЕЙ ОБЛАСТИ ОПРЕДЕЛЕНИЯ ПАРАМЕТРА. Поэтому в компьютерной геометрии чаще всего пользуются не гладкими, а кусочно-гладкими сплайнами, о которых пойдет речь ниже.

ПРОГРАММА 2. Сплайн Лагранжа с возможностью управления опорными точками.

DynamicModule[{pts,pts0,L,t,p,n=6,tt},pts0={{6.0,0.0},{0.0,0.0},{4.8,5.2},{0.3,5.5},{9.0,3.0},{4.3,8.2},{10.0,6.0}};
tt={1.0,2.0,3.0,4.0,5.0,6.0,7.0};
p[i_,s_]:=Apply[Times,(s-tt[[#]])&/@Drop[Range[n+1],{i+1}]];
L=Table[p[i,t]/p[i,tt[[i+1]]],{i,0,n}];
Manipulate[Show[{ParametricPlot[L.pts,{t,tt[[1]],tt[[-1]]},PlotRange->{{-15,15},{-15,15}}],Graphics[{Red,Text[ToString[#-1],pts[[#]]+{1.2,0}]&/@Range[n+1],Text["Сплайн Лагранжа = сплайн Ньютона",{8,-14}]}]}],{{pts,pts0},Locator}]]

.1.2.  СПЛАЙН НЬЮТОНА

СПЛАЙН НЬЮТОНА представляет собой тот же самый сплайн Лагранжа, но в другой форме, удобной для перевычисления в случае добавления новых опорных точек.

Сплайн Ньютона записывается в виде r(t)=(Subscript[a, 0])+(Subscript[a, 1])(t-(Subscript[t, 0]))+(Subscript[a, 2])(t-(Subscript[t, 0]))(t-(Subscript[t, 1]))+…+(Subscript[a, n])(t-(Subscript[t, 0]))…(t-(Subscript[t, n-1])), где векторные коэффициенты (Subscript[a, 0]), ...,(Subscript[a, n]) ищутся из условий прохождения сплайна через точки (Subscript[p, 0]),… ,(Subscript[p, n]). А именно, координаты точек (Subscript[p, 0]),… ,(Subscript[p, n]) определяют коэффициенты (Subscript[a, 0]),… ,(Subscript[a, k]), 0≤ k ≤ n, с помощью системы уравнений:

(Subscript[p, 0])=(Subscript[a, 0]), 
(Subscript[p, 1])=(Subscript[a, 0])+(Subscript[a, 1])((Subscript[t, 1])-(Subscript[t, 0])), 
(Subscript[p, 2])=(Subscript[a, 0])+(Subscript[a, 1])((Subscript[t, 2])-(Subscript[t, 0]))+(Subscript[a, 2])((Subscript[t, 2])-(Subscript[t, 0]))((Subscript[t, 2])-(Subscript[t, 1])), 
.................................................
(Subscript[p, n])=(Subscript[a, 0])+(Subscript[a, 1])((Subscript[t, n])-(Subscript[t, 0]))+(Subscript[a, 2])((Subscript[t, n])-(Subscript[t, 0]))((Subscript[t, n])-(Subscript[t, 1]))+... +(Subscript[a, n])((Subscript[t, n])-(Subscript[t, 0]))...((Subscript[t, n])-(Subscript[t, n-1])).

Эта система уравнений решается за один проход, так как ее матрица является нижне-треугольной. Для любого 1≤ k ≤ n+1 из первых k уравнений определяются первые k неизвестных (Subscript[a, 0 ]),… ,(Subscript[a, k-1]). Если мы добавим еще одну опорную точку (Subscript[p, n+1]), то достаточно будет дополнительно определить еще один коэффициент (Subscript[a, n+1]) из еще одного, добавленного уравнения. Коэффициенты (Subscript[a, 0]) ,… ,(Subscript[a, n]) останутся прежними и их пересчитывать не надо.

ЗАДАЧА. Написать программу для сплайна Ньютона с возможностью добавлять точку.

.1.3.  СПЛАЙНЫ ЭРМИТА

СПЛАЙНЫ ЭРМИТА предназначены для соединения двух заданных кривых новой кривой, полиномиальной по параметру, так, чтобы получившаяся совокупная кривая имела заданную степень гладкости в точках стыковки. Отсюда следует, что сплайн Эрмита задается двумя своими крайними точками и заданным числом m производных радиус-вектора в каждой из крайних точек. Сплайн Эрмита имеет вид полинома степени (2m+1) по t. Степень полинома называется ПОРЯДКОМ СПЛАЙНА. При этом, (2m+2) векторных коэффициентов этого полинома определяются из системы (2m+2) уравнений (m+1 условие для каждой из двух крайних точек).

Сплайн Эрмита может быть определен формулой:

r(t)=(Underoverscript[∑ , (i=0), m])( (Subscript[g, i])(t,m)(Subscript[r, 0]^i)+(Subscript[h, i])(t,m)(Subscript[r, 1]^i) ),

где  (Subscript[r, 0]^( i)), (Subscript[r , 1]^i), 0≤ i≤ m, – заданные значения радиус-вектора кривой и его m производных на концах, а  (Subscript[g, i])(t,m), (Subscript[h, i])(t,m) даются формулами:

(Subscript[g, i])(t,m)    =     (1/(i!))   (∑_(j=0)^(m-i)) (Subscript[C, 2m+1-i ]^j) (t^( i+j)) (((1-t))^( 2m+1-(i+j))),
(Subscript[h, i])(t,m)   =   (((-1))^i)(1/(i!))   (∑_(j=0)^(m-i)) (Subscript[C, 2m+1-i ]^j) (t^( 2m+1-(i+j))) (((1-t))^( i+j)) .

ЗАДАЧА. Написать программу для сплайна Эрмита с возможностью добавлять точку.  Запрограммировав сплайны Эрмита, проверить на Mathematica, что производные радиус вектора равны чему надо.

.1.4.  СОСТАВНОЙ СПЛАЙН ЭРМИТА

СОСТАВНОЙ СПЛАЙН ЭРМИТА предназначен для построения кривой заданной степени гладкости, проходящей через заданную последовательность опорных точек. Он получается путем соединения соседних опорных точек сплайнами Эрмита одного и того же порядка. При этом, значения производных радиус-вектора в опорных точках (требуемые для построения составляющих сплайнов Эрмита) необходимо как-то доопределять. Это можно сделать по-разному, исходя из различных соображений. Простые способы доопределения производных приводят к простым вычислениям, однако при неравномерном расположении опорных точек такой сплайн может образовывать необоснованные петли и изгибы. Более сложные способы доопределения производных (псевдоупругие сплайны Эрмита, например) дают хорошие результаты, но приводят к большим вычислениям.

ПРИМЕР. Составной сплайн Эрмита порядка 3 ( m=1). Для построения такого сплайна необходимо предварительно доопределить первые производные радиус-вектора в опорных точках (Subscript[p, 0 ]),… , (Subscript[p, n]). Для промежуточных и крайних опорных точек способ доопределения различен.

1) Промежуточные опорные точки. В промежуточных опорных точках можно, например, положить производные (Subscript[q, i])=r'((Subscript[t, i])) равными

(Subscript[q, i])=(((Subscript[p, i+1])-(Subscript[p, i-1]))/((Subscript[t, i+1])-(Subscript[t, i-1]))) ,  1≤i≤n-1,

т.е. рассмотреть аппроксимацию производных.

Однако, при неправномерном расположении опорных точек такой способ доопределения производных может привести к появлению необоснованных петель. Поэтому обычно пользуются другими схемами доопределения (Subscript[q, i]). Например, можно положить:

(Subscript[q, i])=(Subscript[s, i+1])(((Subscript[p, i])-(Subscript[p, i-1]))/((Subscript[s, i])+(Subscript[s, i+1]))) +(Subscript[s, i])(((Subscript[p, i+1])-(Subscript[p, i]))/((Subscript[s, i])+(Subscript[s, i+1])))

или

(Subscript[q, i])=(Subscript[s, i])(((Subscript[p, i])-(Subscript[p, i-1]))/((Subscript[s, i])+(Subscript[s, i+1]))) +(Subscript[s, i+1])(((Subscript[p, i+1])-(Subscript[p, i]))/((Subscript[s, i])+(Subscript[s, i+1]))),  1≤i≤n-1,

где (Subscript[s, i])=|(Subscript[p, i])-(Subscript[p, i-1])| – расстояние между соседними опорными точками.

2) Крайние опорные точки. В двух крайних опорных точках (Subscript[p, 0]), (Subscript[p, n]), если в них не требуется гладко стыковать сплайн с другими, уже заданными кривыми, обычно пользуются условием равенства нулю третьих производных радиус-вектора на концах (условие свободных концов кривой):  (r^((3)))((Subscript[t, 0])) = 0,  (r^((3)))((Subscript[t, n])) = 0.

Можно показать, что соответствующие значения первых производных на концах составного сплайна Эрмита третьего порядка даются формулами: (Subscript[q, 0])=2((Subscript[p, 1])-(Subscript[p, 0]))-(Subscript[q, 1]) ,  (Subscript[q, n])=2((Subscript[p, n])-(Subscript[p, n-1]))-(Subscript[q, n-1]).

УПРАЖНЕНИЕ. С помощью программы, рисующей составной сплайн Эрмита третьего порядка, изучить поведение этого сплайна в зависимости от расположения опорных точек, конфигурации узлов и способа доопределения производных в промежуточных точках (поменяйте в программе способ доопределения производных в промежуточных точках на другой и сравните результат).

ПРОГРАММА 3. Составной сплайн Эрмита 3 порядка с 6 управляемыми мышью опорными точками.

DynamicModule[{pts,pts0,n,t,tt,i,j,S,g0,g1,h0,h1,w,q,r},
pts0={{6.0,-5.0},{0.0,0.0},{4.8,5.2},{0.3,7.5},{6.8,12.2},{-2.3,0.0}};
tt={1.0,2.0,3.0,6.0,7.0,8.0};
n=Length[pts0]-1;

Manipulate[
S=Table[{r[i,t,pts],tt[[i]]≤t<=tt[[i+1]]},{i,1,n}];
Show[{ParametricPlot[
Piecewise[S],{t,tt[[1]],tt[[-1]]},PlotRange->{{-10,20},{-10,20}}
],
Graphics[{Red,Text[ToString[#-1],pts[[#]]+{1.2,0}]&/@Range[n+1],Text["Составной сплайн Эрмита 3 порядка",{7,-9}]}]}],{{pts,pts0},Locator}],Initialization:>(
g0[x_]:=2 (x^3)-3 (x^2)+1;
g1[x_]:=(x^3)-2 (x^2)+x;
h0[x_]:=3 (x^2)-2 (x^3);
h1[x_]:=(x^3)-(x^2);
w[i_,t_]:=(t-tt[[i]])/(tt[[i+1]]-tt[[i]]);
q[i_,pts_]:=If[i==1,2(pts[[2]]-pts[[1]])-(pts[[3]]-pts[[1]])/(tt[[3]]-tt[[1]]),If[i==n+1,2(pts[[n+1]]-pts[[n]])-(pts[[n+1]]-pts[[n-1]])/(tt[[n+1]]-tt[[n-1]]),(pts[[i+1]]-pts[[i-1]])/(tt[[i+1]]-tt[[i-1]])]];
r[i_,t_,pts_]:=g0[w[i,t]]*pts[[i]]+g1[w[i,t]]*q[i,pts]+h0[w[i,t]]*pts[[i+1]]+h1[w[i,t]]*q[i+1,pts]
)
]

..  ПСЕВДОУПРУГИЕ СПЛАЙНЫ ЭРМИТА

При построении составного сплайна Эрмита третьего порядка неизвестные производные (Subscript[q, 0]),…  ,(Subscript[q, n]) в точках (Subscript[p, 0]),…  ,(Subscript[p, n]) могут быть определены также из следующих соображений.

Представим себе, что через точки (Subscript[p, 0])=((Subscript[x, 0]),(Subscript[y, 0])),…  ,(Subscript[p, n])=((Subscript[x, n]),(Subscript[y, n])) должна пройти упругая нить, которая стремится занять положение, соответствующее минимуму запасенной энергии ее внутренних сил упругости. Эта энергия, как известно из физики, выражается интегралом от квадрата кривизны k(s) кривой, по которой пролегает нить:

E=(1/2)(∫_( (Subscript[t, 0]))^( (Subscript[t, n ])))c(s)  (k^2)(s)ds ,

где s – натуральный параметр на кривой, а c(s) – коэффициент упругости нити.  Мы будем считать, что c(s) = c = const.

Рассмотрим здесь для простоты случай, когда сплайн может быть задан графиком функции y=y(x). Поскольку минимизация энергии упругих сил является, вообще говоря, достаточно трудоемкой вычислительной задачей, мы заменим ее на более простую, но достаточную для наших целей задачу. А именно, заменим энергию E на ПСЕВДОУПРУГУЮ энергию, заменив в вышеприведенной формуле кривизну k на модуль второй производной |(((d^2)y)/(d (x^2)))|. Это приведит к гораздо более простой задаче минимизации интеграла

(Subscript[E, pseudo])=(1/2) с (∑_(i=0)^(n-1))(Subscript[∫ , (Subscript[x, i])]^(Subscript[x, i+1])) (((y''(x)))^2) dx

за счет выбора величин (Subscript[y, 0]^( ')),… ,(Subscript[y, n]^( ')), где (Subscript[y, i]^( '))=y'((Subscript[x, i])). Условие стационарной точки целевой функции, получаемое приравниванием к нулю ее частных производных по  (Subscript[y, i]^( ')) ,  имеет вид

((∂ E)/(∂ (Subscript[y, j]^( ')))) = с (∑_(i=0)^(n-1))(Subscript[∫ , (Subscript[x, i])]^(Subscript[x, i+1])) y''(x) ((∂ y''(x))/(∂ (Subscript[y, j]^( '))))  dx =0,   j=0,…  ,n.

Это условие может быть приведено к следующему виду

(Subscript[y, i-1, i]^('')) ((Subscript[x, i]))=(Subscript[y, i, i+1]^(''))((Subscript[x, i])), i=1,…  ,n-1,  
(Subscript[y, 0, 1]^('')) ((Subscript[x, 0])) = (Subscript[y, n-1, n]^(''))((Subscript[x, n])) = 0,

где (Subscript[y, i,i+1])(x) – это обозначение y(x) на отрезке [(Subscript[x, i]),(Subscript[x, i+1])].

Другими словами, наше условие означает, что вторые производные y''(x), а следовательно и кривизна кривой, должны быть непрерывны в точках стыковки (Subscript[p, 1]),… ,(Subscript[p, n-1]), а на концах кривой y'' обращается в нуль.

С помощью несложных преобразований можно привести это условие к виду линейной системы уравнений  на (Subscript[y, i]^(')) с ленточной 3-шаговой матрицей:

(Subscript[y, i-1]^('))  (1/(Subscript[Δ, i-1]))+2  (Subscript[y, i]^('))  ((1/(Subscript[Δ, i-1]))+(1/(Subscript[Δ, i])))+  (Subscript[y, i+1]^(')) (1/(Subscript[Δ, i]))  = 3  (((Subscript[y, i])-(Subscript[y, i-1]))/(Subscript[Δ, i-1]^2))+3(((Subscript[y, i+1])-(Subscript[y, i]))/(Subscript[Δ, i]^2)),  1≤i≤n-1;
 2  (Subscript[y, 0]^('))  (Subscript[Δ, 0])+(Subscript[y, 1]^(')) (Subscript[Δ, 0]) = 3((Subscript[y, 1])-(Subscript[y, 0]));     
(Subscript[y, n-1]^('))(Subscript[Δ, n-1])+2(Subscript[y, n]^('))  (Subscript[Δ, n-1]) = 3((Subscript[y, n])-(Subscript[y, n-1])),

где (Subscript[Δ, i])=(Subscript[x, i+1])-(Subscript[x, i]).

Это – система из (n+1)-го уравнения с (n+1) неизвестным (Subscript[y, 0]^(')),… ,(Subscript[y, n]^(')). Поскольку матрица системы имеет ленточный вид, то она допускает численное решение за линейное по n время.

Более подробно с псевдоупругими сплайнами Эрмита можно ознакомится по учебнику Н.Н.Голованова, Д.П.Ильютко, Г.В.Носовского и А.Т.Фоменко "Компьютерная геометрия", разделы 6.1.3 – 6.1.4.

УПРАЖНЕНИЕ. Написать программу пострения псевдоупругого сплайна Эрмита на плоскости, имеющего вид графика функции и проходящего через заданную последовательность опорных точек с возрастающими абсциссами. Сравнить псевдоупругий сплайн Эрмита и составной сплайн Эрмита 3 порядка, проходящий через одни и те же опорные точки путем изображения обоих сплайнов на одном графике.

.1.5.  КУБИЧЕСКИЕ СПЛАЙНЫ

КУБИЧЕСКИЙ СПЛАЙН – это кривая класса (C^2), проходящая через последовательность заданных опорных точек (Subscript[p, 0]),… ,(Subscript[p, n]) при заданных значениях параметра (Subscript[t, 0]),… ,(Subscript[t, n]) и на каждом участке [(Subscript[t, i]),(Subscript[t, i+1])] (то есть между каждой парой соседних опорных точек) имеющая вид кубической кривой.

Это значит, что зависимость ее радиус-вектора от параметра имеет вид многочлена 3 степени на каждом участке, но коэффициенты многочлена могут быть различны на различных участках. Требуется, чтобы кривая в целом имела непрерывные первые и вторые производные радиус-вектора по параметру. Естественно, последнее условие достаточно удовлетворить лишь в промежуточных опорных точках, поскольку во всех других точках кривая является полиномиальной, а следовательно, бесконечно гладкой.

При полиномиальной интерполяции степень полинома в данной постановке задачи не может быть, вообще говоря, сделана меньше 3. Поэтому кубический сплайн определяется однозначно и можно показать, что на каждом отрезке [(Subscript[t, i]),(Subscript[t, i+1])] изменения параметра t он имеет следующий вид:

r(t) = (Subscript[s, i])  (((((Subscript[t, i+1])-t))^3)/(6((Subscript[t, i+1])-(Subscript[t, i]))))+(Subscript[s, i+1]) ((((t-(Subscript[t, i])))^3)/(6((Subscript[t, i+1])-(Subscript[t, i])))) + ( ((Subscript[p, i])/((Subscript[t, i+1])-(Subscript[t, i]))) -(Subscript[s, i])  (((Subscript[t, i+1])-(Subscript[t, i]))/6) )((Subscript[t, i+1])-t)+( ((Subscript[p, i+1])/((Subscript[t, i+1])-(Subscript[t, i]))) -(Subscript[s, i+1]) (((Subscript[t, i+1])-(Subscript[t, i]))/6) )(t-(Subscript[t, i])) .

Векторные коэффициенты (Subscript[s, i]), участвующие в последнем выражении, определяются из условий согласованности первых производных радиус-вектора кривой на концах отрезков [(Subscript[t, i]),(Subscript[t, i+1])], что приводит к следующей системе уравнений для определения (Subscript[s, i]):

(Subscript[s, i-1 ])((Subscript[t, i])-(Subscript[t, i-1]))+2(Subscript[s, i]) ((Subscript[t, i+1])-(Subscript[t, i-1]))+(Subscript[s, i+1 ])((Subscript[t, i+1])-(Subscript[t, i]))  =  6 (((Subscript[p, i+1])-(Subscript[p, i]))/((Subscript[t, i+1])-(Subscript[t, i]))) - 6  (((Subscript[p, i])-(Subscript[p, i-1]))/((Subscript[t, i])-(Subscript[t, i-1]))) ,   i=1,…  ,n-1.

Это система из (n-1) векторногоу равнения (или, что то же самое, 2(n-1) скалярного уравнения) с ленточной трех-диагональной матрицей и (n+1) векторным неизвестным (Subscript[s, 0]), ...,(Subscript[s, n]) (что соответствует 2(n+1) скалярных неизвестных). Если кривая замкнута, то (Subscript[s, 0]) = (Subscript[s, n]) и в систему  добавляется еще одно уравнение для i=n: (Subscript[s, n+1])=(Subscript[s, 1]).

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

Если кривая незамкнута, то необходимо задать граничные условия. Это делается, вообще говоря, неоднозначно. Например, если рассматривать сплайн как гибкую нить, концы которой не имеют нагрузки на изгиб, получим такое условие: (Subscript[s, 0]) = (Subscript[s, n]) = 0.

Еще один способ задания граничных условий следующий: (Subscript[s, 0])=(Subscript[s, 1]),    (Subscript[s, n])=(Subscript[s, n-1 ]).

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

Вышеприведенное выражение для кубического сплайна на участке  [(Subscript[t, i]),(Subscript[t, i+1])] можно переписать и в другом виде, через местный параметр w: r(t)=(1-w)(Subscript[p, i])+w (Subscript[p, i+1])+( (-2w+3(w^2)-(w^3))  (Subscript[s, i])+(-w+(w^3))  (Subscript[s, i+1]) )  (((((Subscript[t, i+1])-(Subscript[t, i])))^2)/6), где w =  (((t - (Subscript[t, i]))/((Subscript[t, i+1]) - (Subscript[t, i]))))  –  местный параметр на участке [(Subscript[t, i]),(Subscript[t, i+1])].

ЗАМЕЧАНИЯ.

В отличии от составного сплайна Эрмита 3 порядка, для кубического сплайна изменение любой из опорных точек всегда приводит к перевычислению всего сплайна.

Для создания сплайнов использются подпакет SplineFit, просто надо загрузить Splines` из Splines Package. Для этого используется функция SplineFit [data, type], которая возвращает сплайн функцию для данных data, используя сплайн-аппроксимацию типа type – по умолчанию это кубический сплайн Cube (другие типы – Bezier и CompositeBezier).

ПРОГРАММА 4. Кубический сплайн с 6-ю управляемыми мышью опорными точками.

DynamicModule[{pts,pts0,n=6,spline},pts0={{6.0,-5.0},{0.0,0.0},{4.8,5.2},{20.3,10.5},{6.8,12.2},{20.0,20.0}};
Manipulate[
Show[{ParametricPlot[spline[pts][x],{x,0,6}],Graphics[{Red,Text[ToString[#-1],pts[[#]]+{1.4,0}]&/@Range[n]}]},PlotRange->{{-1,21},{-6,21}}],{{pts,pts0},Locator}],Initialization:>(Needs["Splines`"];spline[pts_]:=SplineFit[pts,Cubic])]

.2. КРИВЫЕ БЕЗЬЕ

.2.1. ОБЫЧНЫЕ КРИВЫЕ БЕЗЬЕ

Определение. КРИВЫЕ БЕЗЬЕ получаются из следующей постановки задачи. Для любого набора из (n+1) опорных точек (Subscript[p, 0]) ,...,(Subscript[p, n]), требуется
найти такие функции (Subscript[f, 0])(t),...,(Subscript[f, n])(t), не зависящие от (Subscript[p, 0]) ,...,(Subscript[p, n]) и имеющие вид полиномов по t на заданном отрезке изменения параметра (который
мы без ограничения общности можем считать отрезком [0,1]), что кривая r(t)=(Subscript[p, 0])(Subscript[f, 0])(t)+...+(Subscript[p, n])(Subscript[f, n])(t), обладает следующими свойствами:

1) r(0)=(Subscript[p, 0]), r(1)=(Subscript[p, n,])

2) искомая кривая не зависит от выбора декартовой системы координат в (R^n) или, другим словами, вышеприведенная формула, определяющая ее радиус-вектор, инвариантна относительно движений пространства;

3) если все опорные точки (Subscript[p, 0]) ,...,(Subscript[p, n]) лежат в одной плоскости, то искомая кривая – плоская и лежит внутри выпуклой оболочки своих опорных точек (Subscript[p, 0]) ,...,(Subscript[p, n])

Определение. Фиксируем целое n>0.  БАЗИСОМ БЕРНШТЕЙНА на отрезке [0,1] называется система из (n+1) полинома степени n, задаваемых формулами (Subscript[B, i,n])(t)=((Subscript[C, n])^i)(((t^i)(1-t))^(n-i)),0 ≤ i≤ n.

Следующая программа рисует полиномы Бернштейна.

ПРОГРАММА 1. Полиномы Бернштейна.

DynamicModule[{i,k,n,DivZ,PowerZ,Bernstein},
DivZ[s_,0]:=0;
DivZ[s_,t_]:=s/t;
PowerZ[(x_)?((#1 == 0) & ), (y_)?((#1 == 0) & )] := 1; 
PowerZ[x_, y_] := x^y; 
Bernstein[k_, n_, t_]:=Binomial[n,k]*PowerZ[1-t,n-k]*PowerZ[t,k];
Manipulate[Module[{ps,f},If[i>n,i=n];ps=Table[{Hue[.67,.6,.6],Opacity[.1(8+2Boole[i==k])],Thickness[0.001(1+4Boole[k==i])]},{k,0,n}];
f[t_]:=Table[Bernstein[k,n,t],{k,0,n}];Show[Plot[Evaluate@f[x],{x,0,1},PlotStyle->ps,PlotPoints->25,PlotRange->{0,1.1},Filling->If[fill,{i+1->{Axis,{Opacity[.15],Green}}},None],PlotLabel->Style[(Subscript[B, i,n])[t]==Bernstein[i,n,t],14]],Graphics[If[n==0,{},With[{j=i+.3Boole[i==0]-.3Boole[i==n]},Arrow[{{.5,1.1},{DivZ[j,n],Bernstein[i,n,DivZ[j,n]]}}]]]],ImageSize->{450,300},ImagePadding->{{25,25},{25,38}}]],
{{n,12,"Степень полинома n"},0,20,1},{{i,5,"Номер полинома i"},0,n,1},{{fill,True,"Заполнение"},{True,False},ControlPlacement->Bottom}]
]

Оказывается в качестве функций (Subscript[f, 0])(t),...,(Subscript[f, n])(t) можно взять базис Бернштейна.

Определение. Пусть дана (n+1) опорная точка (Subscript[p, 0]) ,...,(Subscript[p, n]). Тогда  КРИВОЙ БЕЗЬЕ, определяемой точками (Subscript[p, 0]) ,...,(Subscript[p, n]), называется следующая кривая (Subscript[r, B])(t)=(∑_(i=0)^n)(Subscript[B, i,n])(t)(Subscript[p, i]). Число n называется СТЕПЕНЬЮ кривой Безье.

ЗАМЕЧАНИЯ.

1) Кривая Безье является полиномиальной, а следовательно бесконечно гладкой на всей области определения параметра [0,1].

2) Построение кривой Безье не требует сложных вычислений. Участвующие в ее определении функциональные коэффициенты (полиномы Бернштейна) зависят лишь от числа опорных точек, но не от их координат. Поэтому кривая Безье легко пересчитывается при любом изменении опорных точек. Такой пересчет требует лишь порядка n арифметических операций.

3) Из определения кривой Безье непосредственно следует, что она проходит через свои крайние опорные точки при значениях параметра t=0 и
t=1 соответственно. Однако через все промежуточые опорные точки кривая Безье, вообще говоря, НЕ ПРОХОДИТ. Таким образом, опорные точки являются лишь СРЕДСТВАМИ УПРАВЛЕНИЯ кривой Безье.

4) В своих крайних опорных точках кривая Безье касается звеньев своей опорной ломаной (ломаной, построенной по опорным точкам).

*** Mathematica имеет встроенную функцию построения кривой Безье BezierCurve[{p1, p2, p3, p4}], однако эта функция строит кривую Безье только
по 4 опорным точкам. Если в качесте ее аргуметов задать больше, чем 4 точки, вмето кривой Безье будет построена составная кривая Безье, имеющая изломы в каждой четвертой опорной точке (попробуйте). Для того, чтобы в Mathematica построить кривую Безье по более, чем 4 опорным точкам, надо пользоваться функцией BezierFunction[pts], где [pts] – список координат опорных точек. Эта функция задает зависимость радиус-вектора кривой Безье от параметра. Кривую Безье можно изобразить, применив к ней оператор  ParametricPlot или ParametricPlot3D.

ЗАМЕЧАНИЕ. Вторая программа демонстрирует пункт 4 замечания.

ПРОГРАММА 2. Кривая Безье, построенная по четырем опорным точкам (точки управляемы).

DynamicModule[{pts,pts0,n=4},
pts0={{-3,0},{0,4},{2,4},{3,2}};
Manipulate[
Show[ {Graphics[{BezierCurve[pts],Green,Line[pts],
Red,Text[ToString[#-1],pts[[#]]+{0.4,0}]&/@Range[4]}]},PlotRange->{{-4,5},{-2,5}}],{{pts,pts0},Locator}]
]

ПРОГРАММА 3. Кривая Безье, построенная по пяти опорным точкам, с управляемыми второй и предпоследней точками. Цель программы – показать, что кривая Безье всегда касается звеньев своей опорной ломанной в крайних точках.

DynamicModule[{u,v,t,pts},
Manipulate[
pts={{-3,0},u,{0,0},v,{3,2}};
BezierF=BezierFunction[pts];
Show[{Graphics[{RGBColor[.12,.34,.65],Line[{{-3,0},u}],Line[{v,{3,2}}],AbsolutePointSize[6],Point[{u,v}],RGBColor[1,.21,0],PointSize[Large],Point[{{-3,0},{0,0},{3,2}}]},Axes-> True],
Graphics[{Red,Text[ToString[#-1],pts[[#]]+{0.3,0}]&/@Range[5],Text["Составной сплайн Эрмита 3 порядка",{7,-9}]}],ParametricPlot[BezierF[t],{t,0,1},PlotStyle-> {{Red,AbsoluteThickness[2]}}]},PlotRange-> 6,ImageSize-> 6 72,Background-> GrayLevel[.98]],{{u,{2,-3}},{-6,-6},{6,6}},{{v,{-6,-3}},{-6,-6},{6,6}},ControlType->Locator]
]

ПРОГРАММА 3. Пример пространственной кривой Безье, построенной по 8 опорным точкам. Все точки управляемы покоординатно с помощью ползунков.

DynamicModule[{pts,p0,p1,p2,p3,p4,p5,p6,p7,t,g},
Manipulate[pts:={p0,p1,p2,p3,p4,p5,p6,p7};
Show[
{
ParametricPlot3D[g[pts][t],{t,0,1},PlotStyle->{Green,Thick},Axes->True,AxesEdge->Automatic,AxesOrigin->{0,0,0}, PlotRange->{{-15,15},{-15,15},{-15,15}}],
Graphics3D[{Red,PointSize[Large],Point[pts[[#]]]&/@Range[8],Text[ToString[#-1],pts[[#]]+{0.5,0.5,0.5}]&/@Range[8]}]
}
],
{{p0,{0.0,0.0,0.0}},{-1,-1,-1},{1,1,1}},{{p1,{2.0,1.0,1.0}},{1,0,0},{3,2,2}},{{p2,{1.8,5.2,2.4}},{0.8,4.2,1.4},{2.8,6.2,3.4}},{{p3,{-2.3,6.5,3.0}},{-3.3,4.5,2.0},{-1.3,7.5,4.0}},{{p4,{-4.0,4.0,5.0}},{-5.0,3.0,4.0},{-3.0,5.0,6.0}},{{p5,{-7.0,1.0,6.0}},{-8.0,0.0,5.0},{-6.0,2.0,7.0}},{{p6,{4.8,2.2,8.4}},{3.8,1.2,7.4},{5.8,3.2,9.4}},{{p7,{9.3,6.5,12.0}},{8.3,5.5,11.0},{10.3,7.5,13.0}},ControlType-> Slider
],Initialization:> (g[pts_]:=BezierFunction[pts]) ]

Для построения кривых Безье существуют разные алгоритм. Следующий пример демонстрирует один из таких алгоритмов – алгоритм де Кастелье, см. Н.Н.Голованов, Д.П.Ильютко, Г.В.Носовский, А.Т.Фоменко, ``Компьютерная геометрия''.

АЛГОРИТМ де КАСТЕЛЬЕ. Пусть даны опорные точки (Subscript[p, 0]) ,...,(Subscript[p, n]). Обозначим (Subscript[r, i])(t;k) кривую Безье, построенную по точкам (Subscript[p, i]) ,...,(Subscript[p, i+k]), i+k≤ n. Тогда искомая кривая Безье – это (Subscript[r, 0])(t;n), и рекуррентная формула де Кастелье для вычисления кривой Безье имеет вид

(Subscript[r, i])(t;k)=(1-t)(Subscript[r, i])(t;k-1)+t (Subscript[r, i+1])(t;k-1), (Subscript[r, i])(t;0)=(Subscript[p, i]).

ПРОГРАММА 4. Алгоритм де Кастелье.

DynamicModule[{s,r,b,t, p},
b[t_,i_,r_]:=(1-t) b[t,i,r-1]+t b[t,i+1,r-1];
b[t_,0,1]=(1-t)p[0]+t p[1];
b[t_,1,1]=(1-t) p[1]+t p[2];
b[t_,2,1]=(1-t) p[2]+t p[3];
p[0]={0.0,0.5};
p[1]={0.5,2.5};
p[2]={3.0,2.0};
p[3]={4.0,0.0};
Manipulate[Module[{plt},plt=ParametricPlot[b[t,0,3],{t,0,s},PlotStyle->{Blue,Thickness[0.01]}][[1]];Which[r==1,Graphics[{{Lighter[Gray],Line[{p[0],p[1],p[2],p[3]}]},{PointSize[0.015],Point[{p[0],p[1],p[2],p[3]}]},{Purple,PointSize[0.015],Point[{b[s,0,1],b[s,1,1],b[s,2,1]}]},{Purple,Thickness[0.005],Line[{p[0],b[s,0,1]}],Line[{p[1],b[s,1,1]}],Line[{p[2],b[s,2,1]}]}},PlotRange->{{-0.1,4.1},{-0.1,2.6}},ImageSize->{300,275}],r==2,Graphics[{{Lighter[Gray],Line[{p[0],p[1],p[2],p[3]}]},{PointSize[0.015],Point[{p[0],p[1],p[2],p[3]}]},{Lighter[Gray],Line[{{b[s,0,1],b[s,1,1]},{b[s,1,1],b[s,2,1]}}]},{Purple,PointSize[0.015],Point[{b[s,0,1],b[s,1,1],b[s,2,1],b[s,0,2],b[s,1,2]}]},{Purple,Thickness[0.005],Line[{{p[0],b[s,0,1]},{p[1],b[s,1,1]},{p[2],b[s,2,1]},{b[s,0,1],b[s,0,2]},{b[s,1,1],b[s,1,2]}}]}},PlotRange->{{-0.1,4.1},{-0.1,2.6}},ImageSize->{300,275}],r==3,Graphics[{plt,{Lighter[Gray],Line[{p[0],p[1],p[2],p[3]}]},{PointSize[0.015],Point[{p[0],p[1],p[2],p[3]}]},{Lighter[Gray],Line[{{b[s,0,1],b[s,1,1]},{b[s,1,1],b[s,2,1]},{b[s,0,2],b[s,1,2]}}]},{Purple,PointSize[0.015],Point[{b[s,0,1],b[s,1,1],b[s,2,1],b[s,0,2],b[s,1,2],b[s,0,3]}]},{Purple,Thickness[0.005],Line[{{p[0],b[s,0,1]},{p[1],b[s,1,1]},{p[2],b[s,2,1]},{b[s,0,1],b[s,0,2]},{b[s,1,1],b[s,1,2]},{b[s,0,2],b[s,0,3]}}]}},PlotRange->{{-0.1,4.1},{-0.1,2.6}},ImageSize->{300,275}]]],{{r,3,"Шаг рекурсии."},{1,2,3}},{{s,0.5,"t (параметр на кривой)"},0.01,1}]
]

ПРОГРАММА 5. Следующий пример строит последовательность кривых Безье в форме ``цветка'', соседние кривые имеют только одну общую крайнюю точку. Все опорные точки расположены на двух окружностях: первая имеет радиус 1, а радиус второй управляется ползунком.

Manipulate[
DynamicModule[{ptlists, img, pts2},
pts2 = 
Table[ If[ EvenQ[i], 1.,r]*{Sin[2 Pi i *iMul/ iDiv], Cos[2 Pi i *iMul/ iDiv]},
{i,0, iLoops/2* LCM[iDiv,iDeg,iMul]}];
ptlists = Partition[ pts2, iDeg+1,iDeg]; 
img = Graphics[{
Table[{If[bColor, Hue[2*k/Length[ptlists]], Black], BezierCurve[ptlists[[k]], SplineDegree->iDeg]},
{k, Length[ptlists]}] },ImageSize->{375,375},Axes->True]],
{{bColor,True,  "Расцветка кривых Безье"}, {True, False}},
{ {iDiv, 23, "Количество равностоящих точек окружностей, в которых (возможно не во всех) будут располагаться все опорные точки"}, 2, 40,1,ImageSize->Tiny,Appearance->"Labeled"},
{{iDeg, 9, "Степень кривых Безье"}, 1, 11,1,ImageSize->Tiny,Appearance->"Labeled"},
{{iMul, 7, "Шаг, с помощью которого выбираются опорные точки"}, 1, 10,1,ImageSize->Tiny,Appearance->"Labeled"},
{{r, .5, "Радиус второй окружности"}, 0, 1.5,.01,ImageSize->Tiny,Appearance->"Labeled"},
{{iLoops,1, "Добавление опорных точек"}, 1, 4,1,ImageSize->Tiny,Appearance->"Labeled"},
ControlPlacement->Left
]

.2.2.  РАЦИОНАЛЬНЫЕ КРИВЫЕ БЕЗЬЕ

Рациональная кривая Безье строится по опорным точкам с заданными ВЕСАМИ. Веса точек – положительные вещественные числа. Рациональная кривая Безье является обобщением обычной кривой Безье и сводится к последней в случае, когда все опорные точки имеют одинаковые веса.

Определение. РАЦИОНАЛЬНОЙ КРИВОЙ БЕЗЬЕ, построенной по опорным точкам (Subscript[p, 0]) ,...,(Subscript[p, n]) с приписанными им положительными вещественными весами (Subscript[ω, 0]),...,(Subscript[ω, n]), называется кривая r(t)=(((∑_(i=0)^n)(Subscript[B, i,n])(t)(Subscript[ω, i])(Subscript[p, i]))/((∑_(i=0)^n)(Subscript[B, i,n])(t)(Subscript[ω, i]))).

В случае равных весов (Subscript[ω, 0])=...=(Subscript[ω, n]) знаменатель в последних формулах превращается в постоянную, а рациональная кривая Безье – в обычную кривую Безье.

ЗАМЕЧАНИЯ.

1) Так же, как и обычная кривая Безье, рациональная кривая Безье всегда находится в выпуклой оболочке своих опорных точек, проходит через первую и последнюю опорные точки и, вообще говоря, не проходит через промежуточные опорные точки.

2) Так же, как и для обычной кривой Безье, в своих крайних опорных точках рациональная кривая Безье касается звеньев своей опорной ломаной (ломаной, построенной по опорным точкам).

3) Так же, как и для обычной кривой Безье, i-тая производная на конце рациональной кривой Безье определяется (i+1) крайней опорной точкой, а также весами этих точек.

4) В отличие от обычных кривых Безье, рациональные кривые Безье способны ТОЧНО (а не приблизительно) представлять дуги кривых второго порядка. Это очень важно для ряда инженерных приложений, в которых требуется, чтобы та или иная кривая в точности являлась дугой, скажем, эллипса, параболы или гиперболы (поскольку от этого зависят физические свойства той детали, в описании которой участвует данная кривая).

5) По сравнению с обычной кривой Безье, рациональная кривая Безье обладает более мощными средствами управления. Ее управляющими параметрами являются теперь уже не только опорные точки, но и их веса. Веса являются принципиально новым множеством управляющих параметров, по сравнению с опорными точками, поскольку  изменение весов и изменение опорных точек действует на рациональную кривую Безье по-разному.

6) При смещении опорной точки (Subscript[p, i]) на вектор (Subscript[ξ, i]) все точки рациональной кривой Безье смещаются в направлении вектора (Subscript[ξ, i]) (на разные величины). Поскольку обычная кривая Безье является частным случаем рациональной кривой Безье (в случае равенства весов), то же верно и для обычной кривой Безье.

7) При увеличение веса (Subscript[ω, i]) опорной точки (Subscript[p, i]) все точки рациональной кривой Безье испытывают смещение в направлении опорной точки (Subscript[p, i]) (на разные величины). Это смещение может оказаться нулевым (например, концы кривой не смещаются при изменении весов), но вектор смещения всегда направлен от исходного положения точки кривой к опорной точке (Subscript[p, i]).

*** Mathematica не имеет встроенной поддержки рациональных кривых Безье. В частности, функция BezierFunction[pts] не работает, если в качестве ее аргумента взять последовательность вещественных чисел (весов), а не координат точек пространства размерности 2 или больше. Поэтому при вычислении знаменателя формулы радиус-вектора рациональной кривой Безье с помощью Mathematica следует дополнить список весов еще одной координатой до двумерных векторов, а затем использовать лишь первую компоненту получившейся двумерной вектор-функции.

ПРОГРАММА 6. Рациональньная кривая Безье, построенная с по 8-ми опорным точкам и их весам, с возможностью непосредственного управления опорными точками с помощью мыши, а весами – с помощью движков, расположенных слева от графика. Для сравнения на том же графике изображена также и обычная кривая Безье, построенная по тем же самым опорным точкам, но без использования весов. Первоначально все веса сделаны одинаковыми, поэтому после запуска программы обе кривые совпадают. При изменении весов они будут расходиться.

DynamicModule[{pts0,w0,w1,w2,w3,w4,w5,w6,w7,ww,n},
pts0={{-10.0,0.0},{-3.0,7.0},{1.1,-7.2},{4.3,-8.5},{-16.0,12.0},{-2.1,1.7},{-8.2,1.4},{12.2,5.4}};
n=Length[pts0];
Manipulate[
ww={w0,w1,w2,w3,w4,w5,w6,w7};
Show[{Graphics[{Red,Text[ToString[#-1],pts[[#]]+{1.2,0.2}]&/@Range[n],Text["Рациональная кривая Безье с весами (Subscript[ω, i]) в сравнении с обычной кривой
     Безье",{-7.2,13}]}],
ParametricPlot[{fBez[pts,ww][t]/fwBez[ww][t][[1]],f0Bez[pts][t]},{t,0,1}]
},
Axes->True,PlotRange->{{-30,15},{-10,15}}],
{{pts,pts0},Locator},
{{w0,1},1,100},{{w1,1},1,100},{{w2,1},1,100},{{w3,1},1,100},{{w4,1},1,100},{{w5,1},1,100},{{w6,1},1,100},{{w7,1},1,100},
Initialization:> (
f0Bez[pts_]:=BezierFunction[pts];
fwBez[ww_]:=BezierFunction[Table[{ww[[i]],1},{i,1,n}]];
fBez[pts_,ww_]:=BezierFunction[pts ww];
)]
]

.2.3. УВЕЛИЧЕНИЕ ЧИСЛА ОПОРНЫХ ТОЧЕК РАЦИОНАЛЬНОЙ КРИВОЙ БЕЗЬЕ БЕЗ ИЗМЕНЕНИЯ КРИВОЙ

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

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

Определим новые опорные точки (Subscript[p', 0]) ,...,(Subscript[p', n+1]) и их веса (Subscript[ω', 0]),...,(Subscript[ω', n+1]) по формулам:

(Subscript[p', i])=((i (Subscript[ω, i-1])(Subscript[p, i-1])+(n+1-i)(Subscript[ω, i])(Subscript[p, i]))/(i (Subscript[ω, i-1])+(n+1-i)(Subscript[ω, i]))),
(Subscript[ω', i])=((i (Subscript[ω, i-1])+(n+1-i)(Subscript[ω, i]))/(n+1)).

Количество новых опорных точек и весов будет на единицу больше, и они будут определять ту же самую рациональную кривую Безье, что и прежде, причем, в той же самой параметризации (см. учебник Н.Н.Голованов, Д.П.Ильютко, Г.В.Носовский, А.Т.Фоменко "Компьютерная геометрия").

Отметим, что неопределенные величины (Subscript[p, n+1]),(Subscript[ω, n+1]),(Subscript[p, -1,])(Subscript[ω, -1]) участвующие в правой части последних формул, входят туда с нулевыми коэффициентами, поэтому доопределять их нет необходимости.

ЗАМЕЧАНИЕ. При добавлении опорной точки весь набор опорных точек (кроме крайних) меняется, а число опорных точек увеличивается на единицу.

ЗАДАЧА. Написать программу для увеличения-уменьшения количества опорных точек рациональной кривой (по 1.2.3.)

.2.4. ДЕЛЕНИЕ КРИВОЙ БЕЗЬЕ НА ДВЕ ЧАСТИ

Для заданного (SuperStar[t]) из интервала (0,1) требуется построить два новых набора опорных точек и весов {(Subscript[p, 0]^a),...,(Subscript[p, n]^a)}, {(Subscript[ω, 0]^a),...,(Subscript[ω, n]^a)} и {(Subscript[p, 0]^b),...,(Subscript[p, n]^b)}, {((Subscript[ω, 0])^b),...,((Subscript[ω, n])^b)} так, чтобы кривые Безье, построенные по этих новым наборам, совпали бы (как геометрическое место точек) соответственно с двумя дугами исходной кривой Безье, на которые она разбивается точкой r((SuperStar[t])).

Данная задача имеет важное значение в теории кривых Безье. Необходимость деления кривых Безье возникает в тех случаях, когда требуется изменить
какую-то часть кривой, оставляя некоторую другую ее часть неизменной. Поскольку смещение любой опорной точки кривой Безье меняет всю ее в целом, возникает необходимость разбивать кривую Безье на несколько частей. Это дает возможность менять только одну из получившихся частей исходной кривой Безье (ставших уже самостоятельными кривыми Безье), не трогая остальные части исходной кривой.

Мы рассмотрим здесь для простоты деление обычной кривой Безье. Деление рациональной кривой Безье может быть сведено к делению обычной кривой Безье (см. учебник "Компьютерная геометрия").

Процедура построения точек {(Subscript[p, 0]^a),...,(Subscript[p, n]^a)} и {(Subscript[p, 0]^b),...,(Subscript[p, n]^b)}, определяющих две новые кривые Безье, делящие исходную на две части, может быть описана следующим образом.
1) Исходный многоугольник с вершинами в опорных точках назовем МНОГОУГОЛЬНИКОМ НУЛЕВОГО УРОВНЯ (Subscript[p, 0]^((0)))=(Subscript[p, 0]),...,(Subscript[p, n]^((0)))=(Subscript[p, n]).

2) Строим многоугольники уровней 1,2,...,n так, чтобы вершины многоугольника следующего уровня делили стороны многоугольника предыдущего
уровня в отношении (SuperStar[t]):(1-(SuperStar[t])): (Subscript[p, i]^((k)))=(1-(SuperStar[t]))(Subscript[p, i]^((k-1)))+(SuperStar[t])(Subscript[p, i+1]^((k-1))), i=0,...,n-k. МНОГОУГОЛЬНИК k УРОВНЯ будет иметь n-k+1 вершин (Subscript[p, 0]^((k))),...,(Subscript[p, n-k]^((k))).

3) Положим: (Subscript[p, i]^a)=(Subscript[p, 0]^((i))).

Аналогично определяются новые опорные точки  {(Subscript[p, 0]^b),...,(Subscript[p, n]^b)} для второй половины деления, только надо брать последовательности последних, а не первых вершин многоугольников k-го уровня.