Подскажите алгоритмы экстраполяции гладких кривых

← →
Unknown user © (2007-06-04 15:22) [0]

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

← →
MBo © (2007-06-04 15:40) [1]

>но она работает только внутри области значений входных данных, то есть экстраполяция невозможна

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

← →
Unknown user © (2007-06-04 15:48) [2]

to MBo

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

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

← →
MBo © (2007-06-04 16:01) [3]

Сплайн на каждом промежутке представляет собой полином (чаще всего кубический). После расчета коэффициентов сплайнов на i-ом промежутке между узлами рассчитывается интерполяционная кривая, т.е. значения полинома P(i) в нужном количестве точек, в диапазоне X(i)..X(i+1). Однако ничего не мешает рассчитать точки этой (кубической) кривой и вне указанного диапазона, хотя чем дальше от него, тем меньше будет похоже на правду. Обычно наибольшая ошибка на краях - при использовании модели свободных концов (вторая производная нулевая), а повысить точность можно, используя предположения о значениях на концах первой производной (наклон) или второй производной.

← →
Unknown user © (2007-06-04 16:08) [4]

я использую данную функцию для рассчета коеф. полиномов. (источник: http://alglib.sources.ru/) как раз с граничным условием нулевой второй производной.

подскажите как рассчитать ctbl за пределами диапазона?

procedure Spline3BuildTable(N : Integer; DiffN : Integer; x : TReal1DArray; y : TReal1DArray; BoundL : Double; BoundR : Double; var ctbl : TReal2DArray); var C : Boolean; E : Integer; G : Integer; Tmp : Double; nxm1 : Integer; I : Integer; J : Integer; DX : Double; DXJ : Double; DYJ : Double; DXJP1 : Double; DYJP1 : Double; DXP : Double; DYP : Double; YPPA : Double; YPPB : Double; PJ : Double; b1 : Double; b2 : Double; b3 : Double; b4 : Double; begin N := N-1; g := (n+1) div 2; repeat i := g; repeat j := i-g; c := True; repeat if x[j]<=x[j+g] then begin c := False; end else begin Tmp := x[j]; x[j] := x[j+g]; x[j+g] := Tmp; Tmp := y[j]; y[j] := y[j+g]; y[j+g] := Tmp; end; j := j-1; until not ((j>=0) and C); i := i+1; until not (i<=n); g := g div 2; until not (g>0); SetLength(ctbl, 4+1, N+1); N := N+1; if DiffN=1 then begin b1 := 1; b2 := 6/(x[1]-x[0])*((y[1]-y[0])/(x[1]-x[0])-BoundL); b3 := 1; b4 := 6/(x[n-1]-x[n-2])*(BoundR-(y[n-1]-y[n-2])/(x[n-1]-x[n-2])); end else begin b1 := 0; b2 := 2*BoundL; b3 := 0; b4 := 2*BoundR; end; nxm1 := n-1; if n>=2 then begin if n>2 then begin dxj := x[1]-x[0]; dyj := y[1]-y[0]; j := 2; while j<=nxm1 do begin dxjp1 := x[j]-x[j-1]; dyjp1 := y[j]-y[j-1]; dxp := dxj+dxjp1; ctbl[1,j-1] := dxjp1/dxp; ctbl[2,j-1] := 1-ctbl[1,j-1]; ctbl[3,j-1] := 6*(dyjp1/dxjp1-dyj/dxj)/dxp; dxj := dxjp1; dyj := dyjp1; j := j+1; end; end; ctbl[1,0] := -b1/2; ctbl[2,0] := b2/2; if n<>2 then begin j := 2; while j<=nxm1 do begin pj := ctbl[2,j-1]*ctbl[1,j-2]+2; ctbl[1,j-1] := -ctbl[1,j-1]/pj; ctbl[2,j-1] := (ctbl[3,j-1]-ctbl[2,j-1]*ctbl[2,J-2])/pj; j := j+1; end; end; yppb := (b4-b3*ctbl[2,nxm1-1])/(b3*ctbl[1,nxm1-1]+2); i := 1; while i<=nxm1 do begin j := n-i; yppa := ctbl[1,j-1]*yppb+ctbl[2,j-1]; dx := x[j]-x[j-1]; ctbl[3,j-1] := (yppb-yppa)/dx/6; ctbl[2,j-1] := yppa/2; ctbl[1,j-1] := (y[j]-y[j-1])/dx-(ctbl[2,j-1]+ctbl[3,j-1]*dx)*dx; yppb := yppa; i := i+1; end; i:=1; while i<=n do begin ctbl[0,i-1] := y[i-1]; ctbl[4,i-1] := x[i-1]; Inc(i); end; end; end;

← →
MBo © (2007-06-04 16:30) [5]

>подскажите как рассчитать ctbl за пределами диапазона
нет, таблицу за пределами диапазона рассчитать нельзя, но после расчета таблицы ты должен использовать функцию вычисления значения сплайна, и полиномы крайних промежутков можно продолжить
http://slil.ru/24463873
цветная кривая слева на участке -1..0 аппроксимирована полиномом, рассчитанным на участке 0..1

SplineCoeff(X, Y, Coef);
Caption := IntToStr((GetTickcount - t) div 100);
for i := 0 to n - 2 do begin
dx := X[i + 1] - X[i];
for j := 0 to 9 do begin
XX := 0.1 * j * dx;
YY := CalcSpline(Y, Coef, xx, i);
Series1.AddXY(X[i] + XX, YY);
if i = 0 then
for j := 0 to 9 do begin
XX := -0.1 * j * dx;
YY := CalcSpline(Y, Coef, xx, i);
Series1.AddXY(X[i] + XX, YY);
end; end;

← →
Unknown user © (2007-06-04 16:45) [6]

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

//TArr -массив параметров обозначающий длину кривой в точках 0..N-1 while Ind2<N do begin D:=Dist(InP[Ind1-1],InP[Ind1]); TArr[Ind2]:=TArr[Ind2-1]+D; Inc(Ind2); end; //построение таблиц коеф. полиномов Spline3BuildTable(N,1,TArr,XArr,0,0,TblX); Spline3BuildTable(N,1,TArr,YArr,0,0,TblY); //экстраполяция (эксперементирую) D:=Dist(InP[0],InP[1]); SetLength(TArr,Length(TArr)+5); //продолжаю на 5 шагов for Ind2:=N to High(TArr) do TArr[Ind2]:=TArr[Ind2-1]+D; T:=0; Ind1:=0; while T<TArr[High(TArr)] do begin P.Status:=0; P.X:=Spline3Interpolate(N,TblX,T); P.Y:=Spline3Interpolate(N,TblY,T); . . . T:=T+Step; end;

← →
MBo © (2007-06-04 16:59) [7]

Ага, теперь вижу - ты строишь двумерную интерполяцию кубическими сплайнами с раздельными сплайнами для X и Y - координат с псевдопараметризацией по номеру точки. Это довольно плохой подход - что-то похожее на правду может получиться при более-менее постоянных длинах дуг в промежутках, но все равно даже для интерполяции нельзя надеяться на адекватное описание формы кривой, т.к. эта процедура предназначена для расчета одномерного сплайна - т.е. функции Y(x)

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

Из простого могу посоветовать сплайны Кэтмулл-Рома или кардинальные сплайны. Они не обеспечивают гладкость второго порядка, но для изолиний это не особо важно.

← →
Иксик © (2007-06-04 20:39) [8]

> Иксик © (04.06.07 20:39) [8]
>
> 2MBo ©
> Хороший вы человек!
>

Кто бы сомневался! Имхо, Борис по числу ценных ответов превосходит тут всех мастеров вместе взятых. А уж по соотношению "полезных"/"потрепательских" ответов - так вообще всех в мире!
:)

to MBo

Спасибо за помощь. После некоторого времени потраченного на эксперименты пришел к выводу, что для экстраполяции допустимо использование пониномов не более 2-го порядка, так как высшие порядки кривых могут вести себя непредсказуемо и давать очень большую погрешность. И совсем неплохие результаты дает обычная линейная экстраполяция :) по крайней мере исключены ошибки, связанные с неправильным изгибом кривой, которые довольно часто случаются при использовании полиномов.

> кубическими сплайнами с раздельными сплайнами для X и Y
> - координат с псевдопараметризацией по номеру точки

параметризация у меня не по номеру точки а по длине кривой в данной точке, отложенной от первой точки. лучшего способа обойти разрывы производной при использовании y=f(x) я не придумал

>параметризация у меня не по номеру точки а по длине кривой в данной точке

Да, действительно, не разглядел, что в TArr. Это неплохой метод.

>неплохие результаты дает обычная линейная экстраполяция
Не пробовал ли параболу построить по двум крайним точкам и углу наклона в крайней точке, согласованному со сплайном?

Подскажите алгоритмы экстраполяции гладких кривых Найти похожие ветки