Как нарисовать кривую Лагранжа через полином?

Question

[EvgenyApMr]

Есть задание: нарисовать кривую catmull-rom с помощью полинома n-го порядка.
На странице есть вариант реализации формулы, алгоритм Эйткена.
Код от источника

for (i = 0; i < N; i++){
     P = 1;
     for (j = 0; j < N; j++)
       if (j != i) P = P*(t-ti[j])/(ti[i] - ti[j]);
   }

Моя реализация на delphi c циклом для рисования, но без итоговой формулы.

while a <> 100 do begin
  Inc( a );
  for i:= 0 to Length( xArray ) do begin
    p:=1;
    for j:= 0 to Length(xArray) do begin
      if j <> i then begin
        p:= p+(a/100-yArray[j])/(yArray[i] - yArray[j]);
      end
    end
  end;
  xt:= //тут должна быть формула для расчёта координаты с применением полинома?
  yt:= //тут должна быть формула для расчёта координаты с применением полинома?
  chart1.SeriesList[1].AddXY(xt, yt);
end;

Как должна выглядеть формула для построения точек?

Чтобы получилась что-то похожее на зеленую линию

Comments

[Модератор]

У вас синтаксические ошибки в коде. Приведите рабочий. Иначе будет считаться заданием с последующим удалением.

Answer 1

[Wataru]

Надо на P домножать координаты i-ой точки.

Comments

[EvgenyApMr]

Можете пояснить?

У меня есть верхний цикл while на 100 итераций. Он отмечает 100 точек для создания плавной кривой. При этом нужно, что бы кривая проходила через точки из массивов xArray : yArray, а их там 5.

То есть будущая формула должна как-то учитывать 5 начальных точек в каждой итерации цикла while

[Константин Цветков]

EvgenyApMr,

как-то учитывать 5 начальных точек

Так укажите как. Или хотя бы нарисуйте.

[Модератор]

EvgenyApMr,

как-то учитывать 5 начальных точек

Это называется "вопрос нужно конкретизировать", иначе удаление.

[Wataru]

EvgenyApMr,
Тогда вам надо xt делать a. При вычислении p надо домножать на yArray[i], а не x. Все p в цикле по i надо будет просуммировать - это и будет значение yt.

[EvgenyApMr]

Константин Цветков, должно получиться что-то похожее на зеленую линию. Она нарисована с помощью сплайнов.

[EvgenyApMr]

Wataru, вы это имели в виду?

While a <> 100 do 
        begin
            a:=a+1;
                 for i:= 0 to Length(xArray) do
                  begin
                  p:=1;
                      for j:= 0 to Length(xArray) do
                      begin
                        if j <> i then
                        begin
                            p:= p+(a/100-yArray[j])/(yArray[i] - yArray[j]);
                        end
                      end
                     end;
                     xt:= a;
                     yt:= p;
                     chart1.SeriesList[1].AddXY(xt, yt);

Рисует только прямую

[Wataru]

EvgenyApMr, Нет. Все разности в p остаются xArray. Но потом надо будет p домножить на yarray[i]. И вот это вот прибавить к yt. Внутри цикла по i. А вот xt вычисляется в цикле по a.

Подумайте сами, вот эта вот формула (t-t[j])(t[i]-t[j]) дает 1 при t=t[i] и 0 для t равных другим индексам. Значит, если ее брать в качестве множителя перед значением точки i, то при x равном xarray[i], именно эта точка возьмется с множителем 1, а все остальные возьмутся с множителем 0.

Т.е. Цикл по a (он же по x), цикл по i, внутри считаем формулу лагранжа для всех j != i. После цикла по j это домножаем на y[i] и суммируем.

[EvgenyApMr]

Wataru, наверно, я вас неправильно понял. Попробовал разные варианты реализации пояснения выше, но рисуются только ровные вертикальные полоски за пределами нужной области. Думаю, проблема с xt.

Может попробуете свой вариант показать в виде кода?

Вот, что получилось у меня ориентируясь на совет

While a <> 100 do
        begin
            yt:=0;
            a:=a+1;
                 for i:= 0 to Length(xArray) do
                  begin
                  p:=1;
                      for j:= 0 to Length(xArray) do
                      begin
                        if j <> i then
                        begin
                            p:= p*(a/100-xArray[j])/(xArray[i] - xArray[j]);
                        end
                      end;
                      yt:= yt+yArray[i]*p;
                     end;
                       xt:= //непонятен порядок расчёта


                 chart1.SeriesList[1].AddXY(xt, yt);

   
        end;

[Wataru]

EvgenyApMr, xt должно быть просто a. Ну и в уравнении, где p домножаете, тоже просто a должно быть, а не a/100.

[EvgenyApMr]

Wataru, примерно так?

While a <> 100 do
        begin
            yt:=0;
            a:=a+1;
                 for i:= 0 to Length(xArray) do
                  begin
                  p:=1;

                      for j:= 0 to Length(xArray) do
                      begin
                        if j <> i then
                        begin
                            p:= p*(a-xArray[j])/(xArray[i] - xArray[j]);
                        end
                      end;
                      yt:= yt+yArray[i]*p;
                     end;
                       xt:= a;

                 chart1.SeriesList[1].AddXY(xt, yt);

Результат выходит более похожим на правду, но явно есть что-то не то

[EvgenyApMr]

Wataru, получается, что если xt = a, то значение не может быть больше 100. На графике последняя точка по x имеет значение 106. При таком ограничении пройти через точку не получится.

[Wataru]

EvgenyApMr, Тогда почему вы работаете только с x до 100, если точка, через которую кривая должна пройти, имеет координату 106?

И вообще, циклы у вас странные. Почему-то от 0, хотя в паскале (это же паскаль?) индексация вроде как с 1 была. Судя по графику, похоже, что последняя точка тупо не учитывается у вас почему-то. А в остальном все правильно - кривая проходит через нужные точки.

[EvgenyApMr]

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

Это delphi.

На всякий случай увеличил а до 105 и 106: кривая уходит дальше, но общая картина остаётся неверной.

Для пяти точек совсем плохо получается:


Нету больше идей, как можно поправить код?

[Wataru]

EvgenyApMr,

но общая картина остаётся неверной

В смысле - не верной? Она верная. Такая полиномиальная кривая обязательно на бесконечности x уйдет в бесконечность по y. Полиномы Лагранжа лишь гарантируют, что кривая пройдет через ключевые точки. А дальше она может идти как угодно.

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

[EvgenyApMr]

Wataru,

За ними продолжать кривую нельзя

Подскажите, получится это реализовать ограничением xt и yt диапазоном между минимальным и максимальным значением в xArray и yArray?

[Wataru]

EvgenyApMr,

Подскажите, получится это реализовать ограничением xt и yt диапазоном между минимальным и максимальным значением в xArray и yArray?

Да.

То есть распределяем 100 точек по нужным координатам, соединяем их и получаем нужную кривую.

Всмысле, у вас в массиве xarray 100 точек? Тогда никакая интерполяция вам не нужна - стройте ломанную по этим точкам. Вот если вы хотите провести через ключевые точки гладкую кривую, то вот тут интерполяция уже нужна. Но тогда вы координаты этой кривой считаете и в других точках, между ключевыми. Чем больше вспомогательных точек (меньше шаг по a), тем точнее и глаже будет график.

[Wataru]

EvgenyApMr, а, нет. Я вас не понял. Вам нужно a менять от минимальной до максимальной координаты в xarray.

[EvgenyApMr]

Wataru, спасибо за помощь! Пока оставлю так.

Смущают огромные колебания. В примерах по ссылке в вопросе всё выглядит намного аккуратнее. Возможно, дело в координатах точек.

[Wataru]

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