Интерполяция полиномами второго порядка графически заданных функций

[HotDogich]

"Интерполяция полиномами второго порядка графически заданных функций." - собственно так звучит тема моей курсовой. Как объяснил преподаватель - точки тыкаем сами - далее рисуется сплайн линия, при том при большом кол-ве точек сплайн должен быть более менее гладкий.
Всё это в Delphi. Делать за меня курсовую не надо. Мне лишь нужна "наводка" по алгоритмам - кои будут между 2мя, 3мя ... N введенных точек прорисовывать между ними сплайн.
Погуглив - понял что существуют разные версии формул/алгоритмов по разным сплайнам, да и по методу тоже какие-то разные.
В общем - дайте/посоветуйте пожалуйста мне алгоритм - который подойдет под поставленную задачу.

[KobolD]

Мне кажется Вы немного путаете интерполяцию и аппроксимацию.
При интерполяции порядок полинома должен совпадать с количеством узлов (точек), а при апроксимации можно задавать степень полинома, но не больше числа узлов.
Собственно интерполяция является частным случаем апроксимации.
Обычно апроксимирующая функция это степенной полином Ai*Xi^i, но можно взять и экспоненциальную Ai*exp(i*Xi).
Общий алгоритм такой:
1 По узлам составляется матрица грамма
2 Методом Гаусса находятся коэффициенты Ai
3 строится апроксимирующая функция с использованиев вычисленных коэфициентов.
Вот тут можно почитать http://matlab.exponenta.ru/spline/book1/10.php

[HotDogich]

Спасибо за ответ, но препод мне объяснил что мне нужно лишь будет строить уравнения типа a*(x^2)+b*x+c=f(x) для каждой линии, условную "гладкость" будет обеспечивать уравнение касательной в конце каждой линии, и коэффициенты с уравнения касательной будут связным звеном между линиями, так есть будет строится слау из
a*(x1^2)+b*x1+c=y1
a*(x2^2)+b*x2+c=y2
2*a*x1+b=y1 уравнение касательной.
пока сам не до конца понял суть, но как то так.
А точнее кажись я до сих пор нифига не понял какой ваще будет алгоритм этого сплайна. Даже книжку одну прочел про интерполяцию, аппроксимацию и сплайны - но ни хера не понятно. =(

[KobolD]

Ну тогда смотри, ты будешь строить параболы для каждой пары точек. Через две точки можно провести бесконечное множество парабол, но если задать условие касательности для первой точки то парабола будет только одна. В твоем примере, для первой точки, касательная равна y1 т.е. она касательна к линии паралельной оси Х.
Решая матрицу
x1^2 x1 1 | y1
x2^2 x2 1 | y2
2*x1 1 0 | K
,где для первых двух точек K=y1
находишь коэффициенты a1 b1 c1
потом строишь матрицу для второй и третей точки
x2^2 x2 1 | y2
x3^2 x3 1 | y3
2*x2 1 0 | K
где К уже равно 2*a1*x2+b1 - потому что касательные для первой и второй параболы во второй точке должна быть равны.
решаешь матрицу находишь коэффициенты a2 b2 c2
строишь матрицу для следующих двух точек где K=2*a2*x3+b2 и т.д.

[HotDogich]

Спасибо большое, только что проверил данный алгоритм, решая слау в моей проге, сделанной уже давно для лабораторок, и найдя коэффициенты в ней, подставил их в прогу из другой лабораторки по решению полиномов - получилось на бумажке нечто вроде сплайна. Осталось только запрограммировать - круто однако что надо лишь прописать чтение тыкнутой точки/точек на канвасе - скопипастить 2 процедурки из лаб - и прорисовать всё это.

[HotDogich]

Дописал прогу, у вас тут ошибка была, или опечатка
"где К уже равно 2*a1*x2+b1" там x1 надо, у меня сначала прога была с x2 - там сплайны росли, бешено росли, с каждый разом все выше/ниже, уходили за края экрана уже после 3-4 точки =)

[HotDogich]

Вот кому интересно =) на форме в Events - онмаусдаун - пишем onclick и далее копипастим код - и вуаля прога должна работать

var
  Form1: TForm1;
  xx: array of integer;
  yy: array of integer;
  ac, bc, cc:extended;
  s: integer = 0; //кол-во тыкнутых точек
implementation

{$R *.dfm}





Procedure Narisovat;
var x1, x2, xn, yx: integer;
begin
x1:=xx[s-1]; x2:=xx[s];
if x2>x1 then begin
  for xn:=x1 to x2 do begin
   yx:=round((ac*xn*xn)+(bc*xn)+cc);
   Form1.Canvas.LineTo(xn,yx);
  end;
end;
if x1>x2 then begin
  for xn:=x1 downto x2 do begin
   yx:=round((ac*xn*xn)+(bc*xn)+cc);
   Form1.Canvas.LineTo(xn,yx);
  end;
end;

end;



Procedure SSolve; //тут будем находить коэффициенты для уравнения параболы
var
x, y, i, f, k, w: integer;
m, x1, x2, y1, y2: real;
e: boolean;
A: array[1..3,1..4] of extended;
begin
x1:=xx[s-1];   x2:=xx[s];
y1:=yy[s-1];   y2:=yy[s];
{составляем СЛАУ            тогда как
 x1^2  x1  1 | y1           A[1,1]  A[1,2]  A[1,3] | A[1,4]
 x2^2  x2  1 | y2           A[2,1]  A[2,2]  A[2,3] | A[2,4]
 2*x1  1   0 | K            A[3,1]  A[3,2]  A[3,3] | A[3,4]
где для первых двух точек K=y1, а далее К уже равно 2*a1*x2+b1}
A[1,1]:=x1*x1;  A[1,2]:=x1;  A[1,3]:=1;  A[1,4]:=y1;
A[2,1]:=x2*x2;  A[2,2]:=x2;  A[2,3]:=1;  A[2,4]:=y2;
A[3,1]:=2*x1;   A[3,2]:=1;   A[3,3]:=0;
if s = 2 then A[3,4]:=y1;
if s > 2 then A[3,4]:=2*ac*x1+bc;

                   repeat
e:=false;
for x:=1 to 3 do begin
  for y:=1 to 4 do begin
    if x=y then begin
//начало исправления предотвращения деления на ноль
      if x<3 then begin
      w:=y;
      f:=0;
        while A[x,y]=0 do begin
          e:=true;
           if w=3 then begin w:=0;
            f:=f+1;
           end else begin
            if w>=0 then begin w:=w+1;
             f:=f+1;
            end;
           end;
           if f>5 then begin
             if A[x,y]=0 then begin
              ShowMessage('имеется нулевой столбец при нахождении коэффициента, хз что дальше делать');
              Application.Terminate;
             end;
           end;
            for i:=1 to 3 do begin
             m:=A[i,y];
             A[i,y]:=A[i,w];
             A[i,w]:=m;
            end;
        end;
      end;
//конец исправления предотвращения деления на ноль
         if x=3 then begin
          if A[x,y]=0 then begin ShowMessage('имеется нулевая строка при нахождении коэффициента, хз что дальше делать');
           Application.Terminate;
          end;
         end;//нулевая последняя строка
           for k:=x+1 to 3 do begin     //ниже главной диагонали
            m:=A[k,y] / A[x,y];
            for i:=1 to 4 do A[k,i]:=A[k,i]-A[x,i]*m;
           end;
       {выше главной диагонали}  if e=false then begin
                                     for k:=1 to x-1 do begin
                                      m:=A[k,y] / A[x,y];
                                      for i:=1 to 4 do A[k,i]:=A[k,i]-A[x,i]*m;
                                     end;
                                 end;
    end;
  end;
end;
           until e=false;
//приводим элементы главной диагонали к единичкам
for x:=1 to 3 do begin
A[x,4]:=A[x,4] / A[x,x];
A[x,x]:=1;
end;
//вывод ответа
ac:=a[1,4];
bc:=a[2,4];
cc:=a[3,4];
Narisovat;
end;








procedure TForm1.onclick(Sender: TObject; Button: TMouseButton;
  Shift: TShiftState; X, Y: Integer);
begin
inc(s);
Form1.canvas.pen.color:=clBlack;
Form1.Canvas.Ellipse(x-3, y-3, x+3, y+3);
setlength(xx,s+1);
setlength(yy,s+1);
xx[s]:=X;
yy[s]:=Y;
if s >= 2 then SSolve;
end;

[KobolD]

Так как у тебя всегда матрица будет 3х3 причем с 1 и нулем на определенных местах то можно один раз на бумажке вывести формулу для нахождения a, b, c и не решать каждый раз матрицу. У меня получилось так (но лучше перепроверить)

a=(y1*(x1-x2)+y2-y1)/(x1+x2)^2
b=y1-2*x1*a
c=y2-b*x2^2

[HotDogich]

По сути да, можно и без всего этого, ведь меньше длина кода - меньше шанс ошибок в нем.
Кстати я уже много чего добавил/изменил в проге.
Только вот программа щас рисует чуток не правильно
например каждая новая тыкнутая точка была правее прошлой - всё норм, одну сделал левее - и например сплайн заканчивался в прошлой точке сверху вниз - тогда новый сплайн не продолжает его плавно - сверху вниз перерастая плавно к новой точке, они сразу снизу вверх продолжает сплайн, ща вот пытаюсь как-нить в коде повернуть эту параболу на 180гр - ток не получается. Пытаюсь повернуть при том - отзеркалив относительно горизонтали и вертикали(сам придумал =))

[KobolD]

чесно говоря не очень понял как там получается на самом деле, но суть проблеммы уловил. Попробуй делать так: если координата Х вновь введеной точки меньше чем координата предыдущей то перед вычислением меняешь точки местами Xi <-> X(i-1),
а в матрице в последнюю строку 2*x 1 0 | K подставляешь значение Х точки которая до перестановки была Х(i-1), т.е. чтобы касательность была к нужной точке.