Форум: "Прочее";
Текущий архив: 2007.07.08;
Скачать: [xml.tar.bz2];
ВнизПодскажите алгоритмы экстраполяции гладких кривых Найти похожие ветки
← →
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]2MBo ©
Хороший вы человек!
← →
Германн © (2007-06-05 01:35) [9]
> Иксик © (04.06.07 20:39) [8]
>
> 2MBo ©
> Хороший вы человек!
>
Кто бы сомневался! Имхо, Борис по числу ценных ответов превосходит тут всех мастеров вместе взятых. А уж по соотношению "полезных"/"потрепательских" ответов - так вообще всех в мире!
:)
← →
Unknown user © (2007-06-07 11:06) [10]to MBo
Спасибо за помощь. После некоторого времени потраченного на эксперименты пришел к выводу, что для экстраполяции допустимо использование пониномов не более 2-го порядка, так как высшие порядки кривых могут вести себя непредсказуемо и давать очень большую погрешность. И совсем неплохие результаты дает обычная линейная экстраполяция :) по крайней мере исключены ошибки, связанные с неправильным изгибом кривой, которые довольно часто случаются при использовании полиномов.
> кубическими сплайнами с раздельными сплайнами для X и Y
> - координат с псевдопараметризацией по номеру точки
параметризация у меня не по номеру точки а по длине кривой в данной точке, отложенной от первой точки. лучшего способа обойти разрывы производной при использовании y=f(x) я не придумал
← →
MBo © (2007-06-07 13:23) [11]>параметризация у меня не по номеру точки а по длине кривой в данной точке
Да, действительно, не разглядел, что в TArr. Это неплохой метод.
>неплохие результаты дает обычная линейная экстраполяция
Не пробовал ли параболу построить по двум крайним точкам и углу наклона в крайней точке, согласованному со сплайном?
Страницы: 1 вся ветка
Форум: "Прочее";
Текущий архив: 2007.07.08;
Скачать: [xml.tar.bz2];
Память: 0.5 MB
Время: 0.046 c