unit StudLib;
type
Point=auto class
public
x,y:real;
end;
type
Spline=class
private
n:integer; // количество узлов интерполяции
aK:array[,] of real; // коэффициенты сплайна
P:array of Point; // узлы интерполяции
procedure MakeSpline;
//
// *** Интерполяция табличной функции кубическим сплайном ***
// Функция Spline возвращает коэффициенты сплайна в массиве.
// Является результатом переработки фортран-программы,
// приведенной в кн. Дж.Форсайт, М.Малькольм, К.Моулер
// "Машинные методы математических вычислений",
// М., "Мир", 1980.
//
begin
if n<2 then Exit;
if n=2 then begin
aK[0,0]:=(P[1].y-P[0].y)/(P[1].x-P[0].x); aK[0,1]:=aK[0,0];
(aK[1,0],aK[1,1],aK[2,0],aK[2,1]):=(0,0,0,0);
Exit
end;
aK[2,0]:=P[1].x-P[0].x; aK[1,1]:=(P[1].y-P[0].y)/aK[2,0];
for var i:=1 to n-2 do begin
aK[2,i]:=P[i+1].x-P[i].x;
aK[0,i]:=2*(aK[2,i-1]+aK[2,i]);
aK[1,i+1]:=(P[i+1].y-P[i].y)/aK[2,i];
aK[1,i]:=aK[1,i+1]-aK[1,i]
end;
aK[0,0]:=-aK[2,0]; aK[0,n-1]:=-aK[2,n-2];
aK[1,0]:=0; aK[1,n-1]:=0;
if n<>3 then begin
aK[1,0]:=aK[1,2]/(P[3].x-P[1].x)-aK[1,1]/(P[2].x-P[0].x);
aK[1,n-1]:=aK[1,n-2]/(P[n-1].x-P[n-3].x)-aK[1,n-3]/(P[n-2].x-P[n-4].x);
aK[1,0]:=aK[1,0]*Sqr(aK[2,0])/(P[3].x-P[0].x);
aK[1,n-1]:=-aK[1,n-1]*Sqr(aK[2,n-2])/(P[n-1].x-P[n-4].x)
end;
for var i:=1 to n-1 do begin
var t:=aK[2,i-1]/aK[0,i-1];
aK[0,i]:=aK[0,i]-t*aK[2,i-1];
aK[1,i]:=aK[1,i]-t*aK[1,i-1]
end;
aK[1,n-1]:=aK[1,n-1]/aK[0,n-1];
for var i:=n-2 downto 0 do
aK[1,i]:=(aK[1,i]-aK[2,i]*aK[1,i+1])/aK[0,i];
aK[0,n-1]:=(P[n-1].y-P[n-2].y)/aK[2,n-2]+
aK[2,n-2]*(aK[1,n-2]+2*aK[1,n-1]);
for var i:=0 to n-2 do begin
aK[0,i]:=(P[i+1].y-P[i].y)/aK[2,i]-
aK[2,i]*(aK[1,i+1]+2*aK[1,i]);
aK[2,i]:=(aK[1,i+1]-aK[1,i])/aK[2,i];
aK[1,i]:=3*aK[1,i]
end;
aK[1,n-1]:=3*aK[1,n-1]; aK[2,n-1]:=aK[2,n-2]
end;
public
constructor (pts:array of Point);
begin
n:=pts.Length;
P:=pts;
SetLength(aK,3,n);
MakeSpline;
end;
function Value(u:real):real;
//
// Вычисляет значение сплайна в заданной точке,
// используя коэффициенты, найденные функцией Spline
// (портированная из фортрана функция Seval)
begin
var i,j:integer;
if u
P[n-1].x then i:=n-2
else begin
j:=n-1;
repeat
var k:=(i+j) div 2;
if u
real; tol:real):real;
//
// *** Нахождение действительного корня функции на интервале изоляции ***
// Считается без проверки, что на интервале [ax;bx] корень есть,
// причем единственный.
// Является результатом переработки одноименной фортран-программы,
// приведенной в кн. Дж.Форсайт, М.Малькольм, К.Моулер
// "Машинные методы математических вычислений",
// М., "Мир", 1980.
//
begin
var Sign:real->real:=x->(x=0?0.0:x<0?-1.0:1.0);
var eps:=1.0;
var tol1:real;
repeat
eps:=eps/2; tol1:=1.0+eps
until tol1<=1.0;
var (fa,fb,a,b):=(f(ax),f(bx),ax,bx);
repeat
var (c,d,fc):=(a,b-a,fa);
var e:=d;
repeat
if Abs(fc)0 then q:=-q;
p:=Abs(p);
if (2*p >= 3*xm*q-Abs(tol1*q)) or (p>=Abs(0.5*e*q)) then (d,e):=(xm,d)
else (e,d):=(d,p/q)
end;
(a,fa):=(b,fb);
if abs(d)>tol1 then b:=b+d
else b+=Abs(tol1)*Sign(xm);
fb:=f(b);
until fb*Sign(fc)>0;
until false
end;
function Quanc8(fun:real->real; a,b,abserr,relerr:real):
(real,real,real,integer);
//
// *** Адаптивная квадратурная программа для численного интегрирования ***
// a, b - границы интервала интегрирования
// abserr - заданная предельная абсолютная ошибка
// relerr - заданная предельная относительная ошибка
// * Возвращаемый кортеж *
// [0] - результат интегрирования (res)
// [1] - оценка абсолютной ошибки результата (errest)
// [2] - индикатор надежности (flag)
// [3] - число точек, в которых вычислялась функция (nofun)
//
// Если flag ненулевой, но мал, результат еще можно принять,
// если велик - программа непригодна для такой функции.
// flag=XXX.YYY, где ХХХ - количество подинтервалов
// длины (b-a)/2^30 с недопустимо большой ошибкой вычисления,
// YYY - часть необработанного основного интервала.
//
// Является результатом переработки одноименной фортран-программы,
// приведенной в кн. Дж.Форсайт, М.Малькольм, К.Моулер
// "Машинные методы математических вычислений",
// М., "Мир", 1980.
//
begin
// этап 1
var (res,errest,flag):=(0.0,0.0,0.0);
var nofun:=0;
var (levmin,levmax,levout,nomax):=(1,30,6,5000);
var nofin:=nomax-8*(levmax-levout+128);
var k:=14175.0;
var (w0,w1,w2,w3,w4):=(3956.0/k,23552.0/k,-3712.0/k,41984.0/k,-18160.0/k);
var (area,cor11):=(0.0,0.0);
if(a=b) then begin Result:=(0.0,0.0,0.0,0); Exit end;
// этап 2
var (lev,nim):=(0,1);
var x:=new real[17];
(x[0],x[16]):=(a,b);
var qprev:=0.0;
var f:=new real[17];
f[0]:=fun(x[0]);
var stone:=(b-a)/16;
x[8]:=(x[0]+x[16])/2;
x[4]:=(x[0]+x[8])/2;
x[12]:=(x[8]+x[16])/2;
x[2]:=(x[0]+x[4])/2;
x[6]:=(x[4]+x[8])/2;
x[10]:=(x[8]+x[12])/2;
x[14]:=(x[12]+x[16])/2;
var j:=2;
while j<=16 do begin f[j]:=fun(x[j]); j+=2 end;
nofun:=9;
// этап 3
var qright:=new real[32];
var fsave:=new real[9,31];
var xsave:=new real[9,31];
while True do begin
x[1]:=(x[0]+x[2])/2; f[1]:=fun(x[1]);
j:=3;
while j<=15 do begin
x[j]:=(x[j-1]+x[j+1])/2; f[j]:=fun(x[j]); j:=j+2
end;
nofun+=8;
var step:=(x[16]-x[0])/16;
var qleft:=(w0*(f[0]+f[8])+w1*(f[1]+f[7])+w2*(f[2]+f[6])+w3*(f[3]+f[5])+
w4*f[4])*step;
qright[lev+1]:=(w0*(f[8]+f[16])+w1*(f[9]+f[15])+w2*(f[10]+f[14])+
w3*(f[11]+f[13])+w4*f[12])*step;
var qnow:=qleft+qright[lev+1];
var qdiff:=qnow-qprev; area+=qdiff;
// Этап 4
var esterr:=Abs(qdiff)/1023;
var tolerr:=Max(abserr,relerr*Abs(area))*(step/stone);
var L70:=false;
if lev>=levmin then begin;
if lev>=levmax then begin
flag+=1; L70:=true
end
else if nofun>nofin then begin
nofin*=2; levmax:=levout; flag+=(b-x[0])/(b-a); L70:=true
end
else if esterr<=tolerr then L70:=true;
end;
// этап 5
if not L70 then begin
nim*=2; lev+=1;
for var i:=1 to 8 do begin
fsave[i,lev]:=f[i+8];
xsave[i,lev]:=x[i+8]
end;
qprev:=qleft;
for var i:=8 downto 1 do begin
f[2*i]:=f[i];
x[2*i]:=x[i]
end;
continue
end;
// этап 7
res+=qnow; errest+=esterr; cor11+=qdiff/1023;
while nim.IsOdd do begin
nim:=nim div 2;
lev-=1
end;
nim+=1;
if lev<=0 then break;
qprev:=qright[lev]; x[0]:=x[16]; f[0]:=f[16];
for var i:=1 to 8 do begin
f[2*i]:=fsave[i,lev];
x[2*i]:=xsave[i,lev]
end;
end;
// этап 8
res+=cor11;
if errest=0 then begin
Result:=(res,errest,flag,nofun); Exit
end;
repeat
var temp:=Abs(res)+errest;
if temp<>Abs(res) then begin
Result:=(res,errest,flag,nofun); Exit
end;;
errest*=2;
until false;
end;
function FMin(ax,bx:real; f:real->real; tol:real):real;
//
// *** Нахождение минимума функции на заданном интервале ***
// tol - это не погрешность вычисленного значения, это скорее
// входное значение параметра, влияющего на длину интервала
// неопределенности, получаемого на выходе. Нет смысла задавать
// значение tol меньше 1e-12, это только замедляет работу.
// Является результатом переработки одноименной фортран-программы,
// приведенной в кн. Дж.Форсайт, М.Малькольм, К.Моулер
// "Машинные методы математических вычислений",
// М., "Мир", 1980.
//
begin
var Sign:(real,real)->real:=(x,y)->Abs(x)*(y>=0?1.0:-1.0);
var c:=0.5*(3-Sqrt(5));
var eps:=1.0;
var tol1:real;
repeat
eps/=2; tol1:=eps+1
until tol1<=1;
eps:=Sqrt(eps);
var (a,b):=(ax,bx);
var v:=a+c*(b-a);
var (w,x,e):=(v,v,0.0);
var fx:=f(x);
var (fv,fw):=(fx,fx);
var xm,u:real;
var d:=0.0;
repeat
xm:=0.5*(a+b);
tol1:=eps*Abs(x)+tol/3;
var tol2:=2*tol1;
if Abs(x-xm)<=(tol2-0.5*(b-a)) then begin Result:=x; Exit end;
if Abs(e)<=tol1 then begin
if x>=xm then e:=a-x else e:=b-x
end
else begin
var r:=(x-w)*(fx-fv);
var q:=(x-v)*(fx-fw);
var p:=(x-v)*q-(x-w)*r;
q:=2*(q-r);
if q>0 then p:=-p;
q:=Abs(q); r:=e; e:=d;
if (Abs(p)>=Abs(0.5*q*r)) or (p<=(q*(a-x))) or (p>=(q*(b-x))) then begin
if x>=xm then e:=a-x else e:=b-x
end
else begin
d:=p/q;
u:=x+d;
if (u-a)-tol1 then u:=x+d;
if Abs(d)fx then begin
if u=x then a:=x else b:=x;
v:=w; fv:=fw; w:=x; fw:=fx; x:=u; fx:=fu
end
until false
end;
function PolRt(xcof:array of real; var ier:integer):array of complex;
//
// Нахождение всех корней полинома с действительными коэффициентами
// методом Ньютона - Рафсона (не выше 36-й степени)
//
// Существенно переработанная для языка PascalABC.NET 3.2
// подпрограмма DPOLRT из библиотеки SSP
// (она же - "Пакет научных программ на языке Фортран-4").
//
// Коэффициенты полинома задаются в порядке возрастания степени.
// ier=0 - ошибок не обнаружено
// ier=1 - мало элементов для построения полинома
// ier=2 - степень полинома превышает 36
// ier=3 - не удалось достигнуть приемлемой точности решения за 500 шагов
// ier=4 - коэффициент при старшей степени полинома нулевой
//
begin
var m:=xcof.Length-2;
if m>35 then begin ier:=2; Exit end;
if xcof[m+1]=0 then begin ier:=4; Exit end;
if m<=0 then begin ier:=1; Exit end;
var SSq:complex->real:=x->Sqr(x.Real)+Sqr(x.Imaginary);
ier:=0;
Result:=new complex[m+1];
var sumsq,temp,alpha:real;
var x:complex;
var L50, L60, L100, L135: boolean;
var (ifit,n2,n,nx,nxx):=(0,0,m,m,m+1);
var kj1:=nxx;
var cof:=xcof.Reverse.ToArray;
while n>=0 do begin
var xo:=cplx(0.00500101,0.01000101);
var ini:=0;
(L50,L60,L135):=(false,false,false);
while true do begin
if not l60 then begin
xo:=cplx(-10*xo.Imaginary,-10*xo.Real); x:=xo; ini+=1
end;
L100:=false;
var xpr:=cplx(0.0,0.0);
for var ict:=1 to 500 do begin
var u:=cplx(cof[n+1],0.0);
var xt:=cplx(1.0,0.0);
var ux:=cplx(0.0,0.0);
if u.Real=0 then begin
(x,L50,L135):=(0.0,true,true); nx-=1; nxx-=1;
break
end;
for var i:=1 to n+1 do begin
temp:=cof[n-i+1];
var xt2:=x*xt;
u+=temp*xt2; ux+=Conjugate(i*temp*xt);
xt:=xt2
end;
sumsq:=SSq(ux);
if sumsq=0.0 then
if ifit=0 then begin L60:=false; continue end
else begin x:=xpr; L50:=true; break end;
var dx:=-u*ux/sumsq; x+=dx;
if Abs(dx.Imaginary)+Abs(dx.Real)<1e-12 then begin
L100:=true; break
end
end;
if not L50 then begin
L100:=L100 or (ifit<>0);
if not L100 then
if ini>=5 then begin ier:=3; Exit end
else begin L60:=false; continue end;
for var l:=0 to nxx do Swap(xcof[kj1-l],cof[l]);
Swap(n,nx);
if ifit=0 then begin
(ifit,xpr,L60):=(1,x,true);
continue
end
else break
end
end;
L50:=false;
if L135 then begin
x:=cplx(x.Real,0.0); (sumsq,alpha,L135):=(0.0,x.Real,false); n-=1;
end
else begin
ifit:=0;
if Abs(x.Imaginary)-1e-10*Abs(x.Real)>=0 then begin
alpha:=2*x.Real; sumsq:=SSq(x); n-=2
end
else begin
x:=cplx(x.Real,0.0); (sumsq,alpha):=(0.0,x.Real); n-=1
end
end;
cof[1]:=cof[1]+alpha*cof[0];
for var l:=1 to n do cof[l+1]+=alpha*cof[l]-sumsq*cof[l-1];
Result[n2]:=x; n2+=1;
if sumsq<>0.0 then begin
Result[n2]:=Conjugate(x); n2+=1
end
end
end;
procedure Solve(a:array[,] of real; b:array of real; ipvt:array of integer);
//
// *** Нахождение решения системы линейных уравнений на основе
// LU-разложения вещественной матрицы посредством функции Decomp.
// Массив b - вектор правых частей систеы.
// Массив ipvt является результатом вызова Decomp.
//
// Является результатом переработки одноименной фортран-программы,
// приведенной в кн. Дж.Форсайт, М.Малькольм, К.Моулер
// "Машинные методы математических вычислений",
// М., "Мир", 1980.
//
begin
var n:=a.RowCount-1;
if n<>0 then begin
for var k:=0 to n-1 do begin
var m:=ipvt[k];
Swap(b[m],b[k]);
for var i:=k+1 to n do b[i]+=a[i,k]*b[k]
end;
for var k:=n downto 1 do begin
b[k]/=a[k,k];
var t:=-b[k];
for var i:=0 to k-1 do begin
b[i]+=a[i,k]*t;
end
end
end;
b[0]/=a[0,0]
end;
function Decomp(a:array[,] of real; var cond:real):array of integer;
//
// *** Вычисление LU-разложения вещественной матрицы посредством
// гауссова исключения и оценка обусловленности матрицы ***
// Используется для решения систем линейных уравнений
// При точной вырожденности cond=MaxReal (примерно 10^308)
// Можно вычислить определитель матрицы, найдя произведение
// диагональных элементов матрицы "а" и умножив его на последний
// элемент массива - результата.
// Для решения системы линейных уравнений надо вызвать процедуру
// Solve, репедав ей вектор правых частей уравнения.
//
// Является результатом переработки одноименной фортран-программы,
// приведенной в кн. Дж.Форсайт, М.Малькольм, К.Моулер
// "Машинные методы математических вычислений",
// М., "Мир", 1980.
//
begin
Result:=new integer[a.ColCount];
var n:=a.RowCount-1;
Result[n]:=1;
if n=0 then begin
if a[0,0]<>0 then cond:=1
else cond:=MaxReal;
Exit
end;
var anorm:=a.Cols.Select(r->r.Select(x->Abs(x)).Sum).Max;
for var k:=0 to n-1 do begin
var m:=k;
for var i:=k+1 to n do
if Abs(a[i,k])>Abs(a[m,k]) then m:=i;
Result[k]:=m;
if m0 then begin
for var i:=k+1 to n do a[i,k]:=-a[i,k]/a[k,k];
for var j:=k+1 to n do begin
Swap(a[m,j],a[k,j]);
if a[k,j]<>0 then
for var i:=k+1 to n do a[i,j]+=a[i,k]*a[k,j]
end
end
end;
var work:=new real[a.ColCount];
for var k:=0 to n do begin
var t:=0.0;
if k<>0 then
for var i:=0 to k-1 do t+=a[i,k]*work[i];
var ek:=1.0;
if t<0.0 then ek:=-1.0;
if a[k,k]=0.0 then begin cond:=MaxReal; Exit end;
work[k]:=-(ek+t)/a[k,k];
end;
for var kb:=0 to n-1 do begin
var k:=n-kb-1;
var t:=0.0;
for var i:=k+1 to n do t+=a[i,k]*work[k];
work[k]:=t;
var m:=Result[k];
if m<>k then Swap(work[m],work[k])
end;
var ynorm:=work.Select(x->Abs(x)).Sum;
Solve(a,work,Result);
var znorm:=work.Select(x->Abs(x)).Sum;
cond:=Max(anorm*znorm/ynorm,1.0)
end;
function Svenn(f:real->real; xn,t:real):(real,real,integer);
//
// Поиск начального интервала неопределенности для локализации
// минимума унимодальной функции методом Свенна.
//
// xn - начальное приближение
// t - шаг
// Результат: кортеж из трех элементов, где первые два представляют
// найденный интервал, если третий элемент нулевой и не представляют
// ничего если начальная точка/шаг выбраны неудачно, либо если функция
// не унимодальна а окрестности заданной точки.
// Примечание: если записать функцию по модулю, можно отыскать интервал
// изоляции корня уравнения f(x)=0
//
begin
var k:=0;
var (x1,x2,x3):=(xn-t,xn,xn+t);
var (y1,y2,y3):=(f(x1),f(x2),f(x3));
var d,a,b:real;
if (y1>=y2) and (y2<=y3) then begin
Result:=(x1,x3,0);
Exit
end
else if (y1<=y2) and (y2>=y3) then begin
Result:=(x1,x3,1);
Exit
end
else if (y1>=y2) and (y2>=y3) then (d,a,x2,k):=(t,x2,x2+t,1)
else (d,b,x2,k):=(-t,x2,x2-t,1);
repeat
x1:=x2+Power(2,k)*d;
(y1,y2):=(f(x1),f(x2));
if (y1=y2;
if d=t then Result:=(a,x1,0)
else Result:=(x1,b,0)
end;
end.