2787 lines
93 KiB
ObjectPascal
2787 lines
93 KiB
ObjectPascal
unit NumLibABC;
|
||
|
||
const &NumLibABCVersion = 'NumLibABC 1.0.2 от 20.03.2024';
|
||
|
||
type
|
||
Point=auto class
|
||
|
||
public
|
||
|
||
x,y:real;
|
||
|
||
end;
|
||
|
||
{$region ApproxCheb}
|
||
type
|
||
/// Аппроксимация табличной функции полиномами Чебышева
|
||
/// по методу наименьших квадратов
|
||
ApproxCheb=class
|
||
|
||
private
|
||
nt:integer; // количество точек
|
||
a:array[,] of real; // расширенная матрица системы (см. MatrCreate)
|
||
|
||
procedure MatrCreate;
|
||
// заполняет матрицу а для решения СЛАУ с целью получения
|
||
// коэффициентов интерполяционного полинома степени n
|
||
begin
|
||
a:=new real[r+1,r+2]; // инициализация нулями
|
||
var p,s:real;
|
||
for var i:=0 to r do begin
|
||
(p,s):=(1.0,x[i]); a[i,r+1]:=f[i];
|
||
for var j:=0 to r do begin a[i,j]:=p; p*=s end
|
||
end
|
||
end;
|
||
|
||
procedure Gauss;
|
||
// Базируется на программе 3.1Р
|
||
// В кн.: Мудров А.Е. Численные методы для ПЭВЫМ на языках
|
||
// Бейсик, Фортран и Паскаль. - Томск: МП "Раско", 1991.
|
||
begin
|
||
var n1:=r+1;
|
||
c:=new real[n1];
|
||
var p,s:real;
|
||
for var k:=0 to r do begin
|
||
var k1:=k+1; s:=a[k,k]; var jj:=k;
|
||
for var i:=k1 to r do begin
|
||
p:=a[i,k];
|
||
if Abs(r)>Abs(s) then (s,jj):=(p,i)
|
||
end;
|
||
if jj<>k then
|
||
for var i:=k to n1 do Swap(a[k,i],a[jj,i]);
|
||
for var j:=k1 to n1 do a[k,j]/=s;
|
||
for var i:=k1 to r do begin
|
||
p:=a[i,k];
|
||
for var j:=k1 to n1 do a[i,j]-=a[k,j]*p
|
||
end
|
||
end;
|
||
for var i:=r downto 0 do begin
|
||
s:=a[i,n1];
|
||
for var j:=i+1 to r do s-=a[i,j]*c[j];
|
||
c[i]:=s
|
||
end
|
||
end;
|
||
|
||
public
|
||
|
||
x:array of real; // значения аргументов
|
||
y:array of real; // значения функции
|
||
tol:real; // погрешносмть аппроксимации
|
||
r:integer; // степень аппроксимирующего полинома (r<nt)
|
||
f:array of real; // аппроксимированные значения функции
|
||
c:array of real; // коэффициенты аппроксимирующего полинома
|
||
|
||
constructor (x,y:array of real; eps:real);
|
||
begin
|
||
nt:=x.Length;
|
||
Assert(nt=y.Length,
|
||
'Approx: несовпадение количества точек в аргументе и функции');
|
||
Self.x:=x;
|
||
Self.y:=y;
|
||
f:=ArrFill(nt,0.0);
|
||
Cheb(eps)
|
||
end;
|
||
|
||
/// Аппроксимация табличной функции полиномами Чебышева (МНК)
|
||
procedure Cheb(e:real);
|
||
begin
|
||
r:=0;
|
||
var n:=nt-1;
|
||
var f1:=ArrFill(nt,0.0);
|
||
var f2:=ArrFill(nt,1.0);
|
||
var eps:=ArrFill(nt,0.0);
|
||
while true do begin
|
||
var sum1:=y.Zip(f2,(p,q)->p*q).Sum;
|
||
var sum2:=f2.Select(p->p*p).Sum;
|
||
var al:=sum1/sum2;
|
||
f:=f.Zip(f2,(p,q)->p+al*q).ToArray;
|
||
for var i:=0 to n do eps[i]:=Abs(y[i]-f[i]);
|
||
tol:=Sqrt((eps.Select(t->t*t).Sum)/n);
|
||
if (tol<e) or (r>=n+1) then break;
|
||
sum1:=x.Zip(f2,(p,q)->p*q*q).Sum;
|
||
sum2:=f2.Select(p->p*p).Sum;
|
||
var be:=-(sum1/sum2);
|
||
var g:real;
|
||
if r<>0 then begin
|
||
(sum1,sum2):=(0.0,0.0);
|
||
for var i:=0 to n do begin sum1+=x[i]*f1[i]*f2[i]; sum2+=Sqr(f1[i]) end;
|
||
g:=-(sum1/sum2);
|
||
end;
|
||
r+=1;
|
||
for var i:=0 to n do begin
|
||
var a:=f1[i]; f1[i]:=f2[i]; f2[i]:=(x[i]+be)*f1[i]+g*a
|
||
end
|
||
end
|
||
end;
|
||
|
||
/// Получение коэффициентов найденного полинома в его степенном представлении
|
||
procedure MakeCoef;
|
||
begin
|
||
MatrCreate;
|
||
Gauss
|
||
end;
|
||
|
||
end;
|
||
{$endregion}
|
||
|
||
{$region Decomp}
|
||
type
|
||
/// LU-разложение вещественной матрицы системы линейных уравнений
|
||
Decomp=class
|
||
|
||
private
|
||
|
||
n:integer; // число строк/столбцов исходной матрицы минус 1
|
||
|
||
procedure Decomp;
|
||
//
|
||
// *** Вычисление LU-разложения вещественной матрицы посредством
|
||
// гауссова исключения и оценка обусловленности матрицы ***
|
||
// Используется для решения систем линейных уравнений
|
||
// При точной вырожденности cond=MaxReal (примерно 10^308)
|
||
// Можно вычислить определитель матрицы, найдя произведение
|
||
// диагональных элементов матрицы "а" и умножив его на последний
|
||
// элемент массива - результата.
|
||
// Для решения системы линейных уравнений надо вызвать процедуру
|
||
// Solve, передав ей вектор правых частей уравнения.
|
||
//
|
||
// Является результатом переработки одноименной фортран-программы,
|
||
// приведенной в кн. Дж.Форсайт, М.Малькольм, К.Моулер
|
||
// "Машинные методы математических вычислений",
|
||
// М., "Мир", 1980.
|
||
//
|
||
begin
|
||
Assert(a.RowCount=a.ColCount,'DECOMP: матрица не квадратная');
|
||
ipvt[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;
|
||
ipvt[k]:=m;
|
||
if m<>k then ipvt[n]:=-ipvt[n];
|
||
Swap(a[m,k],a[k,k]);
|
||
if a[k,k]<>0 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[n+1];
|
||
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:=ipvt[k];
|
||
if m<>k then Swap(work[m],work[k])
|
||
end;
|
||
var ynorm:=work.Select(x->Abs(x)).Sum;
|
||
Self.Solve(work);
|
||
var znorm:=work.Select(x->Abs(x)).Sum;
|
||
cond:=Max(anorm*znorm/ynorm,1.0);
|
||
det:=real(ipvt.Last);
|
||
for var i:=0 to a.RowCount-1 do det*=a[i,i];
|
||
end;
|
||
|
||
public
|
||
|
||
a:array[,] of real; // исходная матрица
|
||
cond:real;
|
||
ipvt:array of integer; // вспомогательный вектор
|
||
det:real; // детерминант матрицы a
|
||
|
||
constructor(a:array[,] of real);
|
||
begin
|
||
n:=a.RowCount-1;
|
||
Assert((n=a.ColCount-1) and (n>=0),'DECOMP: исходная матрица не квадратная');
|
||
Self.a:=new real[n+1,n+1];
|
||
for var i:=0 to n do
|
||
for var j:=0 to n do
|
||
Self.a[i,j]:=a[i,j];
|
||
Self.ipvt:=new integer[n+1];
|
||
Decomp
|
||
end;
|
||
|
||
/// Решение системы линейных уравнений для заданного
|
||
/// вектора правых частей уравнения.
|
||
procedure Solve(b:array of real);
|
||
//
|
||
// Массив b - вектор правых частей систеы.
|
||
// Массив ipvt является результатом вызова Decomp.
|
||
//
|
||
// Является результатом переработки одноименной фортран-программы,
|
||
// приведенной в кн. Дж.Форсайт, М.Малькольм, К.Моулер
|
||
// "Машинные методы математических вычислений",
|
||
// М., "Мир", 1980.
|
||
//
|
||
begin
|
||
Assert(b.Length-1=n,'DECOMP/SOLVE: неверная длина вектора правых частей');
|
||
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;
|
||
b[0]/=a[0,0]
|
||
end;
|
||
|
||
end;
|
||
{$endregion}
|
||
|
||
{$region Factors}
|
||
type
|
||
/// Разложение полинома с целочисленными коэффициентами
|
||
/// на рациональные линейные множители вида ux-v
|
||
Factors=class
|
||
|
||
private
|
||
a:array of integer;
|
||
n:integer;
|
||
|
||
public
|
||
|
||
constructor(params a:array of integer);
|
||
begin
|
||
Self.a:=a;
|
||
Self.n:=a.Length-1;
|
||
end;
|
||
|
||
/// Разложение полинома с целочисленными коэффициентами
|
||
/// на рациональные линейные множители вида ux-v
|
||
function Factorize:array[,] of integer;
|
||
//
|
||
// Коэффициенты полинома должны задаваться в порядке возрастания степени.
|
||
// Возвращаемый массив:
|
||
// [0,0] - количество найденных линейных множителей;
|
||
// [0,1] - максимальный НОД коэффициентов полинома;
|
||
// Каждая последующая i-я строка содержит коэффициенты множителя
|
||
// [i,0]=u, [i,1]=v
|
||
//
|
||
// Замечание: В случае, если свободный член нулевой, предварительно следует
|
||
// вынести аргумент за скобку и производить поиск в оставшемся полиноме.
|
||
//
|
||
// Является переработкой алгоритма 75б (язык Алгол-60)
|
||
// В кн. Агеев М.И., Алик В.П., Марков Ю.И.
|
||
// Библиотека алгоритмов 51б-100б. (Справочное пособие).
|
||
// Вып.2. М., "Сов.радио", 1976
|
||
//
|
||
begin
|
||
var f,g,q:integer;
|
||
var (r,c):=(0,1);
|
||
Result:=new integer[n+1,2];
|
||
var b:=a.Reverse.ToArray;
|
||
while b[n]=1 do begin // исключим множители (1*х-0)
|
||
n-=1; r+=1;
|
||
(Result[r,0],Result[r,1]):=(1,0);
|
||
end;
|
||
var Fin:=false;
|
||
for var p:=1 to Abs(b[0]) do
|
||
if b[0] mod p=0 then
|
||
for var q1:=1 to Abs(b[n]) do begin
|
||
q:=q1;
|
||
if q<>1 then begin
|
||
var flag:=true;
|
||
for var j:=0 to n do
|
||
if b[j] mod q<>0 then begin flag:=false; break end;
|
||
if flag then
|
||
for var j:=0 to n do b[j]:=b[j] div q
|
||
end;
|
||
repeat
|
||
(f,g):=(b[0],1);
|
||
for var i:=1 to n do begin g*=p; f:=f*q+g*b[i] end;
|
||
if f=0 then begin
|
||
r+=1;
|
||
(Result[r,0],Result[r,1]):=(p,q);
|
||
for var i:=0 to n do begin f:=(b[i]+f) div p; b[i]:=f; f*=q end;
|
||
n-=1;
|
||
if n<>0 then continue;
|
||
Fin:=true
|
||
end;
|
||
if not Fin then q:=-q;
|
||
until (q>0) or Fin;
|
||
if Fin then break;
|
||
end;
|
||
if n=0 then c*=b[0];
|
||
(Result[0,0],Result[0,1]):=(r,c);
|
||
SetLength(Result,r+1,2)
|
||
end;
|
||
|
||
end;
|
||
{$endregion}
|
||
|
||
{$region FMin}
|
||
type
|
||
/// Нахождение минимума функции на заданном интервале
|
||
FMin=class
|
||
|
||
private
|
||
|
||
a,b,t:real;
|
||
f:real->real;
|
||
|
||
procedure FMin;
|
||
|
||
// Нахождение локального минимума функции f(x) на интервале [a;b]
|
||
// Метод использует комбинацию поиска золотого сечения и последовательной
|
||
// параболической интерполяции. Сходимость никода не бывает хуже, чем при
|
||
// фибоначчиевом поиске. Если функция имеет непрерывную положительную вторую
|
||
// производную в точне минимума, не совпадающем с текущими границами интервала
|
||
// поиска, сходимость сверхлинейная и обычно имеет порядок 1.324...
|
||
// Бесполезно пытаться задавать точность е меньше значения квадратного корня
|
||
// из машинной точности, принимаемого по умолчанию.
|
||
// Результат: кортеж двух переменных типа real:
|
||
// найденные значения аргумента и вычисленной для него функции.
|
||
//
|
||
// На основе версии Fortran90 (John Burkardt).
|
||
// Алгоритм:
|
||
// Richard Brent,
|
||
// Algorithms for Minimization Without Derivatives,
|
||
// Dover, 2002,
|
||
// ISBN: 0-486-41998-3,
|
||
// LC: QA402.5.B74.
|
||
|
||
begin
|
||
var d,fu,u,tol1:real;
|
||
var c:=0.5*(3-Sqrt(5.0)); // квадрат величины, обратной золотому сечению
|
||
var(sa,sb):=(a,b);
|
||
x:=sa+c*(b-a);
|
||
var w:=x;
|
||
var (v,e,fx):=(w,0.0,f(x));
|
||
var fw:=fx;
|
||
var fv:=fw;
|
||
var eps:=1.0;
|
||
repeat
|
||
eps:=0.5*eps;
|
||
tol1:=1.0+eps
|
||
until tol1<=1.0;
|
||
eps:=Sqrt(eps);
|
||
repeat
|
||
var m:=0.5*(sa+sb);
|
||
var tol:=eps*Abs(x)+t;
|
||
var t2:=2*tol;
|
||
// проверка критерия выхода
|
||
if Abs(x-m) <= t2-0.5*(sb-sa) then break;
|
||
// приемлема ли парабола?
|
||
var r:=0.0;
|
||
var q:=r;
|
||
var p:=q;
|
||
if tol<Abs(e) then begin
|
||
r:=(x-w)*(fx-fv); q:=(x-v)*(fx-fw);
|
||
p:=(x-v)*q-(x-w)*r; q:=2.0*(q-r);
|
||
if q>0.0 then p:=-p;
|
||
q:=Abs(q); r:=e; e:=d
|
||
end;
|
||
if (Abs(p)<Abs(0.5*q*r)) and (q*(sa-x)<p) and (p<q*(sb-x)) then begin
|
||
// шаг параболической интерполяции
|
||
d:=p/q; u:=x+d;
|
||
// f не следует вычислять слишком близко к a или b.
|
||
if ((u-sa)<t2) or ((sb-u)<t2) then d:=x<m?tol:-tol;
|
||
end
|
||
else begin
|
||
// шаг золотого сечения
|
||
e:=x<m?sb-x:sa-x; d:=c*e
|
||
end;
|
||
// f не следует вычислять слишком близко к найденному х
|
||
if tol<=Abs(d) then u:=x+d
|
||
else if d>0.0 then u:=x+tol
|
||
else u:=x-tol;
|
||
fu:=f(u);
|
||
// Обновить найденные параметры
|
||
if fu<=fx then begin
|
||
if u<x then sb:=x else sa:=x;
|
||
(v,fv):=(w,fw);
|
||
(w,fw):=(x,fx);
|
||
(x,fx):=(u,fu)
|
||
end
|
||
else begin
|
||
if u<x then sa:=u else sb:=u;
|
||
if (fu<=fw) or (w=x) then begin
|
||
(v,fv):=(w,fw);
|
||
(w,fw):=(u,fu)
|
||
end
|
||
else
|
||
if (fu<=fv) or (v=x) or (v=w) then (v,fv):=(u,fu)
|
||
end
|
||
until false;
|
||
Value:=fx
|
||
end;
|
||
|
||
public
|
||
|
||
x:real; // значение аргумента
|
||
Value:real; // значение функции
|
||
|
||
constructor(f:real->real; a,b:real; t:real:=1.05e-8);
|
||
begin
|
||
Self.a:=a;
|
||
Self.b:=b;
|
||
Self.t:=t;
|
||
Self.f:=f;
|
||
FMin
|
||
end;
|
||
|
||
end;
|
||
{$endregion}
|
||
|
||
{$region FMinN}
|
||
type
|
||
/// минимизация функции многих переменных
|
||
FMinN=class
|
||
|
||
private
|
||
|
||
delta:array of real; // инициализируется HJ
|
||
fbefore:real; // значение функции при входе
|
||
newx:array of real; // инициализируется HJ
|
||
|
||
fun:function(x:array of real):real; // целевая функция
|
||
|
||
function BestNearby:real;
|
||
// Поиск точки, лучше исходной (по одному параметру за шаг)
|
||
// За основу взят текст программы BEST_NEARBY.
|
||
// Источник:
|
||
// TOMS178
|
||
// Minimization by Hooke-Jeeves Direct Search
|
||
// FORTRAN-90 version by John Burkardt.
|
||
// ( Dept. of Scientific Computing, Florida State University ),
|
||
// на условиях GNU LGPL
|
||
begin
|
||
Result:=fbefore;
|
||
var z:=Copy(newx);
|
||
for var i:=0 to n-1 do begin
|
||
z[i]:=newx[i]+delta[i];
|
||
var ftmp:=fun(z);
|
||
if ftmp<Result then Result:=ftmp
|
||
else begin
|
||
delta[i]:=-delta[i];
|
||
z[i]:=newx[i]+delta[i];
|
||
ftmp:=fun(z);
|
||
if ftmp<Result then Result:=ftmp
|
||
else z[i]:=newx[i]
|
||
end
|
||
end;
|
||
newx:=Copy(z);
|
||
end;
|
||
|
||
function IsEqual(x1,x2:array of real; eps:real):Boolean;
|
||
begin
|
||
Result:=true;
|
||
for var i:=0 to x1.Length-1 do
|
||
Result:=Result and (Abs(x1[i]-x2[i])<=eps);
|
||
end;
|
||
|
||
public
|
||
|
||
n:integer;
|
||
iter:integer; // число итераций
|
||
x:array of real; // вектор аргументов
|
||
|
||
constructor (px:array of real; f:function(x:array of real):real);
|
||
begin
|
||
n:=px.Length;
|
||
x:=Copy(px);
|
||
fun:=f;
|
||
iter:=0;
|
||
Randomize;
|
||
end;
|
||
|
||
/// Нахождение безусловного минимума методом адаптивного случайного поиска
|
||
procedure ARS(R:real; var t:real; MaxIter:integer:=100);
|
||
//
|
||
// t - радиус области неопределенности найденных координат;
|
||
// R - минимальная величина шага по аргументам, R<t;
|
||
// MaxIter - максимальное количество итераций
|
||
//
|
||
// Алгоритм приведен в книге Пантелеев А.В., Летова Т.А.
|
||
// Методы оптимизации в примерах и задачах. - М., Высш.шк., 2005.
|
||
//
|
||
begin
|
||
var (a,b,n):=(1.618,0.618,x.Length);
|
||
var (y,z,ksi):=(new real[n],new real[n],new real[n]);
|
||
var (k,j,M):=(0,1,3*n);
|
||
Randomize;
|
||
t:=1.0;
|
||
while true do begin
|
||
var ka:=0.0;
|
||
for var i:=0 to n-1 do begin
|
||
ksi[i]:=2*Random-1.0;
|
||
ka+=Sqr(ksi[i])
|
||
end;
|
||
ka:=Sqrt(ka);
|
||
for var i:=0 to n-1 do y[i]:=x[i]+t*ksi[i]/ka;
|
||
var fx:=fun(x);
|
||
if fun(y)<fx then begin // шаг удачный
|
||
for var i:=0 to n-1 do z[i]:=x[i]+a*(y[i]-x[i]);
|
||
if fun(z)<fun(x) then begin // направление поиска удачное
|
||
x:=Copy(z); t:=a*t; k+=1;
|
||
if k<MaxIter then begin j:=1; continue end
|
||
else break
|
||
end
|
||
end;
|
||
if j<M then j+=1
|
||
else
|
||
if t>R then begin t:=b*t; j:=1 end
|
||
else break
|
||
end
|
||
end;
|
||
|
||
/// Поиск минимума функции методом Монте-Карло
|
||
procedure MKSearch(a,b:array of real; // векторы границ
|
||
var y:real; // лучшее значение функции fun
|
||
m:integer:=1000); // число случайных точек
|
||
// Авторская реализация. Алгоритм тривиален.
|
||
// Найденный вектор параметров сохраняется в свойстве класса x
|
||
// Значение функции для этого вектора возвращается в параметре y
|
||
begin
|
||
var n:=x.Length;
|
||
var tx:=new real[n];
|
||
var ty:real;
|
||
loop m do begin
|
||
for var i:=0 to n-1 do tx[i]:=Random*(b[i]-a[i])+a[i];
|
||
ty:=fun(tx);
|
||
if ty<y then begin x:=Copy(tx); y:=ty end
|
||
end
|
||
end;
|
||
|
||
/// Алгоритм прямого поиска Хука - Дживса
|
||
function HJ(eps:real:=1e-6;
|
||
rho:real:=0.5;
|
||
itermax:integer:=5000):array of real;
|
||
|
||
// За основу взят текст программы HOOKE.
|
||
// Источник
|
||
// TOMS178
|
||
// Minimization by Hooke-Jeeves Direct Search
|
||
// FORTRAN-90 version by John Burkardt.
|
||
// ( Dept. of Scientific Computing, Florida State University ),
|
||
// на условиях GNU LGPL
|
||
// Утверждается что aлгоритм Хука - Дживса хорошо работает
|
||
// при rho=0.5, но плохо при rho=0.6 и лучше при rho=0.8
|
||
|
||
begin
|
||
newx:=Copy(x);
|
||
var xbefore:=Copy(x);
|
||
delta:=x.Select(t->t=0?rho:rho*Abs(t)).ToArray;
|
||
var steplength:=rho;
|
||
iter:=0;
|
||
fbefore:=fun(newx);
|
||
var newf:=fbefore;
|
||
while (iter<itermax) and (eps<steplength) do begin
|
||
iter+=1;
|
||
// поиск новой "лучшей точки" по одной координате за шаг
|
||
newx:=Copy(xbefore);
|
||
newf:=BestNearby;
|
||
// если улучшение получено, идем в этом направлении
|
||
var keep:=1;
|
||
while (newf<fbefore) and (keep=1) do begin
|
||
for var i:=0 to n-1 do begin
|
||
if newx[i]<=xbefore[i] then delta[i]:=-Abs(delta[i])
|
||
else delta[i]:=Abs(delta[i]);
|
||
var tmp:=xbefore[i];
|
||
xbefore[i]:=newx[i];
|
||
newx[i]+=newx[i]-tmp
|
||
end;
|
||
fbefore:=newf;
|
||
newf:=BestNearby;
|
||
if fbefore<=newf then break; // если улучшить не получается...
|
||
keep:=0;
|
||
for var i:=0 to n-1 do
|
||
if 0.5*Abs(delta[i])<Abs(newx[i]-xbefore[i]) then begin
|
||
keep:=1; break
|
||
end;
|
||
end;
|
||
if (eps<=steplength) and (fbefore<=newf) then begin
|
||
steplength*=rho;
|
||
delta.Transform(x->x*rho);
|
||
end
|
||
end;
|
||
Result:=Copy(xbefore);
|
||
end;
|
||
|
||
/// Алгоритм наилучшей пробы с направляющим гиперквадратом
|
||
procedure BPHS(a,b:array of real; // векторы границ
|
||
var y:real; // лучшее значение
|
||
k:integer:=100; // число случайных точек
|
||
m:integer:=1000); // число случайных точек для MKSearch
|
||
//
|
||
// Использует MKSearch. Функция вычисляется k*m раз.
|
||
//
|
||
// Рейзлин В.И. Численные методы оптимизации: учебное пособие
|
||
// / В.И. Рейзлин; Томский политехнический университет.
|
||
// -Томск: Изд-во Томского политехнического университета, 2011.
|
||
|
||
begin
|
||
var alpha:=2.0;
|
||
y:=real.MaxValue;
|
||
var ty:=y;
|
||
var ta:=Copy(a);
|
||
var tb:=Copy(b);
|
||
var brk:=false;
|
||
loop k do begin
|
||
MKSearch(ta,tb,ty,m);
|
||
if ty<y then begin
|
||
y:=ty;
|
||
for var i:=0 to n-1 do begin
|
||
a[i]:=ta[i]; b[i]:=tb[i];
|
||
var h:=(tb[i]-ta[i])/2/alpha;
|
||
brk:=brk or (h+1=1);
|
||
ta[i]:=x[i]-h; brk:=brk or (h+1=1); tb[i]:=x[i]+h;
|
||
end;
|
||
if brk then break
|
||
end
|
||
end
|
||
end;
|
||
|
||
/// Выбор результатов по лучшим из p проб, сделанных BPHS
|
||
function BestP(a,b:array of real;
|
||
eps:real:=0.01; //- версия 1.0.1
|
||
p:integer:=10;
|
||
k:integer:=100;
|
||
m:integer:=1000):List<(real,array of real)>;
|
||
|
||
// Авторская программа. Использует BPHS, MKSearch (косвенно)
|
||
|
||
begin
|
||
var L1:=new List<(real,array of real)>;
|
||
var L2:=new List<(real,array of real)>;
|
||
Result:=new List<(real,array of real)>;
|
||
var a1:=ArrFill(p,true);
|
||
if p>0 then begin
|
||
var i:=0;
|
||
loop p do begin
|
||
var y:=real.MaxValue;
|
||
var ta:=Copy(a);
|
||
var tb:=Copy(b);
|
||
BPHS(ta,tb,y,k,m);
|
||
L1.Add((y,x));
|
||
i+=1;
|
||
end;
|
||
end;
|
||
for var i:=0 to L1.Count-2 do begin
|
||
if a1[i] then begin
|
||
L2.Clear;
|
||
L2.Add(L1[i]);
|
||
for var j:=i+1 to L1.Count-1 do
|
||
if a1[j] then
|
||
if IsEqual(L1[i][1],L1[j][1],eps) then begin
|
||
L2.Add(L1[j]);
|
||
a1[j]:=false
|
||
end;
|
||
var q:=L2.Select(x->x[0]).ToArray.IndexMin;
|
||
if L2.Count>1 then Result.Add(L2[q]) // от случайностей
|
||
end
|
||
end
|
||
end;
|
||
|
||
end;
|
||
{$endregion}
|
||
|
||
{$region Fraction}
|
||
type
|
||
/// Класс обыкновенных дробей произвольной точности
|
||
Fraction=class
|
||
|
||
private
|
||
|
||
class function ToReal(m:BigInteger):(real,integer);
|
||
begin
|
||
var s:=m.ToString;
|
||
var r:real;
|
||
if Real.TryParse(s,r) then Result:=(r,0)
|
||
else Result:=(s.Left(16).ToReal,s.Length-16)
|
||
end;
|
||
|
||
/// НОД пары чисел
|
||
class function GCD(a, b: BigInteger): BigInteger;
|
||
begin
|
||
if a < 0 then
|
||
a := -a;
|
||
while b <> 0 do
|
||
(a, b) := (b, a mod b);
|
||
Result := a
|
||
end;
|
||
|
||
public
|
||
numerator:BigInteger;
|
||
denominator:BigInteger;
|
||
|
||
constructor (num:BigInteger:=0; denom:BigInteger:=1);
|
||
begin
|
||
Assert(denom>0,'Fraction: знаменатель дроби должен быть положительным');
|
||
var n:=GCD(num,denom);
|
||
numerator:=num div n;
|
||
denominator:=denom div n
|
||
end;
|
||
|
||
class function operator +(a,b:fraction):fraction;
|
||
begin
|
||
var n:=GCD(a.denominator,b.denominator);
|
||
Result:=new fraction(a.numerator*(b.denominator div n)+
|
||
b.numerator*(a.denominator div n),a.denominator div n*b.denominator);
|
||
end;
|
||
|
||
class procedure operator +=(var r:fraction; a:fraction);
|
||
begin
|
||
r:=r+a;
|
||
end;
|
||
|
||
class function operator +(a:BigInteger; b:fraction):fraction;
|
||
begin
|
||
Result:=new fraction(a*b.denominator+b.numerator,b.denominator);
|
||
var n:=GCD(Result.numerator,Result.denominator);
|
||
if n<>1 then begin
|
||
Result.numerator:=Result.numerator div n;
|
||
Result.denominator:=Result.denominator div n
|
||
end
|
||
end;
|
||
|
||
class function operator +(a:fraction; b:BigInteger):fraction;
|
||
begin
|
||
Result:=new fraction(b*a.denominator+a.numerator,a.denominator);
|
||
var n:=GCD(Result.numerator,Result.denominator);
|
||
if n<>1 then begin
|
||
Result.numerator:=Result.numerator div n;
|
||
Result.denominator:=Result.denominator div n
|
||
end
|
||
end;
|
||
|
||
class procedure operator +=(var r:fraction; a:BigInteger);
|
||
begin
|
||
r:=r+a;
|
||
end;
|
||
|
||
class function operator -(a,b:fraction):fraction;
|
||
begin
|
||
var n:=GCD(a.denominator,b.denominator);
|
||
Result:=new fraction(a.numerator*(b.denominator div n)-
|
||
b.numerator*(a.denominator div n),a.denominator div n*b.denominator);
|
||
end;
|
||
|
||
class procedure operator -=(var R:fraction; a:fraction);
|
||
begin
|
||
r:=r-a;
|
||
end;
|
||
|
||
class function operator -(a:BigInteger; b:fraction):fraction;
|
||
begin
|
||
Result:=new fraction(a*b.denominator-b.numerator,b.denominator);
|
||
var n:=GCD(Result.numerator,Result.denominator);
|
||
if n<>1 then begin
|
||
Result.numerator:=Result.numerator div n;
|
||
Result.denominator:=Result.denominator div n
|
||
end
|
||
end;
|
||
|
||
class function operator -(a:fraction; b:BigInteger):fraction;
|
||
begin
|
||
Result:=new fraction(a.numerator-b*a.denominator,a.denominator);
|
||
var n:=GCD(Result.numerator,Result.denominator);
|
||
if n<>1 then begin
|
||
Result.numerator:=Result.numerator div n;
|
||
Result.denominator:=Result.denominator div n
|
||
end
|
||
end;
|
||
|
||
class procedure operator -=(var r:fraction; a:BigInteger);
|
||
begin
|
||
r:=r-a;
|
||
end;
|
||
|
||
class function operator *(a,b:fraction):fraction;
|
||
begin
|
||
Result:=new fraction(a.numerator*b.numerator,a.denominator*b.denominator);
|
||
var n:=GCD(Result.numerator,Result.denominator);
|
||
if n<>1 then begin
|
||
Result.numerator:=Result.numerator div n;
|
||
Result.denominator:=Result.denominator div n
|
||
end
|
||
end;
|
||
|
||
class procedure operator *=(var r:fraction; a:fraction);
|
||
begin
|
||
r:=r*a;
|
||
end;
|
||
|
||
class function operator *(a:BigInteger; b:fraction):fraction;
|
||
begin
|
||
Result:=new fraction(a*b.numerator,b.denominator);
|
||
var n:=GCD(Result.numerator,Result.denominator);
|
||
if n<>1 then begin
|
||
Result.numerator:=Result.numerator div n;
|
||
Result.denominator:=Result.denominator div n
|
||
end
|
||
end;
|
||
|
||
class function operator *(a:fraction; b:BigInteger):fraction;
|
||
begin
|
||
Result:=new fraction(a.numerator*b,a.denominator);
|
||
var n:=GCD(Result.numerator,Result.denominator);
|
||
if n<>1 then begin
|
||
Result.numerator:=Result.numerator div n;
|
||
Result.denominator:=Result.denominator div n
|
||
end
|
||
end;
|
||
|
||
class procedure operator *=(var r:fraction; a:BigInteger);
|
||
begin
|
||
r:=r*a;
|
||
end;
|
||
|
||
class function operator /(a,b:fraction):fraction;
|
||
begin
|
||
Result:=new fraction(a.numerator*b.denominator,a.denominator*b.numerator);
|
||
var n:=GCD(Result.numerator,Result.denominator);
|
||
if n<>1 then begin
|
||
Result.numerator:=Result.numerator div n;
|
||
Result.denominator:=Result.denominator div n
|
||
end
|
||
end;
|
||
|
||
class procedure operator /=(var r:fraction; a:fraction);
|
||
begin
|
||
r:=r/a;
|
||
end;
|
||
|
||
class function operator /(a:BigInteger; b:fraction):fraction;
|
||
begin
|
||
Result:=new fraction(a*b.denominator,b.numerator);
|
||
var n:=GCD(Result.numerator,Result.denominator);
|
||
if n<>1 then begin
|
||
Result.numerator:=Result.numerator div n;
|
||
Result.denominator:=Result.denominator div n
|
||
end
|
||
end;
|
||
|
||
class function operator /(a:fraction; b:BigInteger):fraction;
|
||
begin
|
||
Result:=new fraction(a.numerator,a.denominator*b);
|
||
var n:=GCD(Result.numerator,Result.denominator);
|
||
if n<>1 then begin
|
||
Result.numerator:=Result.numerator div n;
|
||
Result.denominator:=Result.denominator div n
|
||
end
|
||
end;
|
||
|
||
class procedure operator /=(var r:fraction; a:BigInteger);
|
||
begin
|
||
r:=r/a;
|
||
end;
|
||
|
||
class function operator -(a:fraction):=
|
||
new fraction(-a.numerator,a.denominator);
|
||
|
||
class function operator =(a:fraction; b:fraction):=(a-b).numerator=0;
|
||
|
||
class function operator =(a:BigInteger; b:fraction):=(a-b).numerator=0;
|
||
|
||
class function operator =(a:fraction; b:BigInteger):=(a-b).numerator=0;
|
||
|
||
class function operator >(a:fraction; b:fraction):=(a-b).numerator>0;
|
||
|
||
class function operator >(a:BigInteger; b:fraction):=(a-b).numerator>0;
|
||
|
||
class function operator >(a:fraction; b:BigInteger):=(a-b).numerator>0;
|
||
|
||
class function operator >=(a:fraction; b:fraction):=(a-b).numerator>=0;
|
||
|
||
class function operator >=(a:BigInteger; b:fraction):=(a-b).numerator>=0;
|
||
|
||
class function operator >=(a:fraction; b:BigInteger):=(a-b).numerator>=0;
|
||
|
||
class function operator <(a:fraction; b:fraction):=(a-b).numerator<0;
|
||
|
||
class function operator <(a:BigInteger; b:fraction):=(a-b).numerator<0;
|
||
|
||
class function operator <(a:fraction; b:BigInteger):=(a-b).numerator<0;
|
||
|
||
class function operator <=(a:fraction; b:fraction):=(a-b).numerator<=0;
|
||
|
||
class function operator <=(a:BigInteger; b:fraction):=(a-b).numerator<=0;
|
||
|
||
class function operator <=(a:fraction; b:BigInteger):=(a-b).numerator<=0;
|
||
|
||
class function operator <>(a:fraction; b:fraction):=(a-b).numerator<>0;
|
||
|
||
class function operator <>(a:BigInteger; b:fraction):=(a-b).numerator<>0;
|
||
|
||
class function operator <>(a:fraction; b:BigInteger):=(a-b).numerator<>0;
|
||
|
||
/// Обратная величина дроби
|
||
function Inv:fraction;
|
||
begin
|
||
Result:=new fraction;
|
||
Result.numerator:=denominator;
|
||
Result.denominator:=numerator
|
||
end;
|
||
|
||
/// Абсолютная величина дроби
|
||
function Abs:fraction;
|
||
begin
|
||
Result:=new fraction;
|
||
Result.numerator:=PABCSystem.Abs(numerator);
|
||
Result.denominator:=denominator
|
||
end;
|
||
|
||
function ToReal:real;
|
||
begin
|
||
var a:=Fraction.ToReal(numerator);
|
||
var b:=Fraction.ToReal(denominator);
|
||
Result:=a[0]/b[0]*Power(10,a[1]-b[1])
|
||
end;
|
||
|
||
function Print(s:string:=' '):fraction;
|
||
begin
|
||
Result:=Self;
|
||
if denominator=1 then Write(numerator,s)
|
||
else Write(numerator,'/',denominator,s)
|
||
end;
|
||
|
||
function ToString:string; override;
|
||
begin
|
||
Result:=numerator.ToString;
|
||
if denominator<>1 then Result+='/'+denominator.ToString
|
||
end;
|
||
|
||
end;
|
||
|
||
/// Обертка для конструктора класса
|
||
function Frc(a:BigInteger):fraction;
|
||
begin
|
||
Result:=new fraction(a,1);
|
||
end;
|
||
|
||
/// Обертка для конструктора класса
|
||
function Frc(a,b:BigInteger):fraction;
|
||
begin
|
||
Result:=new fraction(a,b);
|
||
end;
|
||
|
||
/// Обертка для конструктора класса
|
||
function Frc(a,b,c:BigInteger):fraction;
|
||
begin
|
||
Result:=new fraction(b,c);
|
||
Result:=a+Result
|
||
end;
|
||
{$endregion}
|
||
|
||
{$region Vector}
|
||
type
|
||
/// В Е К Т О Р Ы
|
||
Vector=class
|
||
|
||
public
|
||
|
||
Value:array of real;
|
||
Length:integer;
|
||
|
||
constructor (params px:array of real);
|
||
begin
|
||
Value:=PABCSystem.Copy(px);
|
||
Length:=Value.Length;
|
||
Assert(Value.Length<>0,'Vector: нулевая длина запрещена')
|
||
end;
|
||
|
||
constructor (pa,pb:array of real);
|
||
begin
|
||
Value:=PABCSystem.Copy(pb);
|
||
Length:=Value.Length;
|
||
for var i:=0 to Length-1 do Value[i]-=pa[i];
|
||
Assert(Value.Length<>0,'Vector: нулевая длина запрещена');
|
||
end;
|
||
|
||
constructor (pn:integer);
|
||
begin
|
||
Value:=new real[pn];
|
||
Length:=pn
|
||
end;
|
||
|
||
class function operator + (va,vb:Vector):Vector;
|
||
begin
|
||
Assert(va.Length=vb.Length,'Vector:длины несовместимы');
|
||
Result:=new Vector(va.Value.Zip(vb.Value,(p,q)->p+q).ToArray);
|
||
end;
|
||
|
||
class function operator - (va,vb:Vector):Vector;
|
||
begin
|
||
Assert(va.Length=vb.Length,'Vector:длины несовместимы');
|
||
Result:=new Vector(va.Value.Zip(vb.Value,(p,q)->p-q).ToArray);
|
||
end;
|
||
|
||
class function operator - (va:Vector):Vector; // унарный минус
|
||
begin
|
||
Result:=new Vector(va.Value);
|
||
for var i:=0 to Result.Length-1 do Result.Value[i]:=-Result.Value[i];
|
||
end;
|
||
|
||
class function operator * (va,vb:Vector):real; // скалярное произведение
|
||
begin
|
||
Assert(va.Length=vb.Length,'Vector:длины несовместимы');
|
||
Result:=va.Value.Zip(vb.Value,(p,q)->p*q).Sum;
|
||
end;
|
||
|
||
class function operator * (ra:real; vb:Vector):Vector;
|
||
begin
|
||
Result:=new Vector(vb.Value);
|
||
for var i:=0 to Result.Length-1 do Result.Value[i]*=ra
|
||
end;
|
||
|
||
class function operator * (va:Vector; rb:real):=rb*va;
|
||
|
||
class function operator = (va,vb:Vector):boolean;
|
||
begin
|
||
Assert(va.Length=vb.Length,'Vector:длины несовместимы');
|
||
Result:=(va.Value[0]/vb.Value[0]=1.0);
|
||
var i:=1;
|
||
while Result and (i<=va.Length-1) do begin
|
||
Result:=Result and (va.Value[i]/vb.Value[i]=1.0);
|
||
i+=1
|
||
end
|
||
end;
|
||
|
||
/// векторное произведение
|
||
function VP(vb:Vector):Vector;
|
||
begin
|
||
Result:=new Vector(3);
|
||
Assert(Length=vb.Length,'Vector.VP: длины векторов несовместимы');
|
||
if (Length=2) or (Length=3) then begin
|
||
var r:=new real[3,3] ((1,1,1),(0,0,0),(0,0,0));
|
||
for var j:=0 to Length-1 do begin
|
||
r[1,j]:=Value[j];
|
||
r[2,j]:=vb.Value[j]
|
||
end;
|
||
Result.Value[0]:=r[1,1]*r[2,2]-r[1,2]*r[2,1];
|
||
Result.Value[1]:=-r[1,0]*r[2,2]+r[1,2]*r[2,0];
|
||
Result.Value[2]:=r[1,0]*r[2,1]-r[1,1]*r[2,0]
|
||
end
|
||
else begin
|
||
Assert(false,'Vector.VP: Векторы должны быть двух или трехмерными');
|
||
Result.Value:=ArrFill(3,MaxReal)
|
||
end
|
||
end;
|
||
|
||
/// смешанное произведение
|
||
function MP(vb,vc:Vector):real;
|
||
begin
|
||
Assert((Length=vb.Length) or (Length=vc.Length),
|
||
'Vector.MP: длины векторов несовместимы');
|
||
if (Length=2) or (Length=3) then begin
|
||
var r:=new real[3,3];
|
||
for var j:=0 to Length-1 do begin
|
||
r[0,j]:=Value[j];
|
||
r[1,j]:=vb.Value[j];
|
||
r[2,j]:=vc.Value[j]
|
||
end;
|
||
Result:=r[0,0]*(r[1,1]*r[2,2]-r[1,2]*r[2,1])-
|
||
r[0,1]*(r[1,0]*r[2,2]-r[1,2]*r[2,0])+
|
||
r[0,2]*(r[1,0]*r[2,1]-r[1,1]*r[2,0])
|
||
end
|
||
else begin
|
||
Assert(false,'Vector.VP: Векторы должны быть двух или трехмерными');
|
||
Result:=MaxReal
|
||
end
|
||
end;
|
||
|
||
/// модуль вектора
|
||
function ModV:real;
|
||
begin
|
||
Result:=Sqrt(Value.Select(t->t*t).Sum);
|
||
end;
|
||
|
||
/// орт вектора
|
||
function Ort:Vector;
|
||
begin
|
||
Result:=new Vector(Length);
|
||
var invd:=1.0/Sqrt(Value.Select(t->t*t).Sum);
|
||
for var i:=0 to Length-1 do Result.Value[i]:=Value[i]*invd
|
||
end;
|
||
|
||
/// коллинеарность векторов
|
||
function IsCollinear(vb:Vector):boolean;
|
||
begin
|
||
Assert(Length=vb.Length,'Vector.IsCollinear: длины векторов несовместимы');
|
||
var k:=Value[0]/vb.Value[0];
|
||
Result:=true;
|
||
var i:=1;
|
||
while Result and (i<=Length-1) do begin
|
||
Result:=Result and (Value[i]/vb.Value[i]=k);
|
||
i+=1
|
||
end
|
||
end;
|
||
|
||
/// компланарность векторов
|
||
function IsCoplanar(vb,vc:Vector):boolean;
|
||
begin
|
||
Assert((Length=vb.Length) or (Length=vc.Length),
|
||
'Vector.IsCoplanar: длины векторов несовместимы');
|
||
Result:=Self.MP(vb,vc)=0.0
|
||
end;
|
||
|
||
/// возвращает копию вектора
|
||
function Copy:Vector;
|
||
begin
|
||
Result:=new Vector(Value)
|
||
end;
|
||
|
||
function Print(delim:string:=' '):Vector;
|
||
begin
|
||
Value.Print(delim);
|
||
Result:=new Vector(Value)
|
||
end;
|
||
|
||
function Println(delim:string:=' '):Vector;
|
||
begin
|
||
Value.Println(delim);
|
||
Result:=new Vector(Value)
|
||
end;
|
||
|
||
end;
|
||
{$endregion}
|
||
|
||
{$region Matrix}
|
||
type
|
||
/// М А Т Р И Ц Ы
|
||
Matrix=class
|
||
|
||
public
|
||
|
||
Value: array[,] of real;
|
||
RowCount: integer;
|
||
ColCount: integer;
|
||
|
||
constructor (pa:array[,] of real);
|
||
begin
|
||
(RowCount,ColCount):=(pa.RowCount,pa.ColCount);
|
||
Value:=new real[RowCount,ColCount];
|
||
for var i:=0 to RowCount-1 do
|
||
for var j:=0 to ColCount-1 do
|
||
Value[i,j]:=pa[i,j]
|
||
end;
|
||
|
||
constructor (nRows,nCols:integer);
|
||
begin
|
||
Value:=new real[nRows,nCols];
|
||
RowCount:=nRows;
|
||
ColCount:=nCols
|
||
end;
|
||
|
||
constructor (nRows,nCols:integer; params pa:array of real);
|
||
begin
|
||
Assert(pa.Length=nRows*nCols,
|
||
'Создание MarixR: Неверное количество параметров');
|
||
(RowCount,ColCount):=(nRows,nCols);
|
||
Value:=new real[RowCount,ColCount];
|
||
for var i:=0 to RowCount-1 do
|
||
for var j:=0 to ColCount-1 do
|
||
Value[i,j]:=pa[ColCount*i+j]
|
||
end;
|
||
|
||
class function operator + (pa,pb:Matrix):Matrix;
|
||
// сложение матриц
|
||
begin
|
||
Assert((pa.RowCount=pb.RowCount) and (pa.ColCount=pb.ColCount),
|
||
'Сложение MarixR: Несовпадение размеров');
|
||
Result:=new Matrix(pa.RowCount,pa.ColCount);
|
||
for var i:=0 to pa.RowCount-1 do
|
||
for var j:=0 to pa.ColCount-1 do
|
||
Result.Value[i,j]:=pa.Value[i,j]+pb.Value[i,j]
|
||
end;
|
||
|
||
class function operator - (pa:Matrix):Matrix;
|
||
// смена знака матрицы (унарный минус)
|
||
begin
|
||
Result:=new Matrix(pa.RowCount,pa.ColCount);
|
||
for var i:=0 to pa.RowCount-1 do
|
||
for var j:=0 to pa.ColCount-1 do
|
||
Result.Value[i,j]:=-pa.Value[i,j]
|
||
end;
|
||
|
||
class function operator - (pa,pb:Matrix):Matrix;
|
||
// вычитание матриц
|
||
begin
|
||
Assert((pa.RowCount=pb.RowCount) and (pa.ColCount=pb.ColCount),
|
||
'Вычитание MarixR: Несовпадение размеров');
|
||
Result:=new Matrix(pa.RowCount,pa.ColCount);
|
||
for var i:=0 to pa.RowCount-1 do
|
||
for var j:=0 to pa.ColCount-1 do
|
||
Result.Value[i,j]:=pa.Value[i,j]-pb.Value[i,j]
|
||
end;
|
||
|
||
class function operator * (r:real; pb:Matrix):Matrix;
|
||
// умножение скаляра на матрицу
|
||
begin
|
||
Result:=new Matrix(pb.RowCount,pb.ColCount);
|
||
for var i:=0 to pb.RowCount-1 do
|
||
for var j:=0 to pb.ColCount-1 do
|
||
Result.Value[i,j]:=r*pb.Value[i,j]
|
||
end;
|
||
|
||
class function operator * (pa:Matrix; r:real):=r*pa;
|
||
// умножение матрицы на скаляр
|
||
|
||
class function operator * (pa,pb:Matrix):Matrix;
|
||
// умножение матриц
|
||
begin
|
||
Assert(pa.ColCount=pb.RowCount,
|
||
'Matrix: умножение матриц. Несовпадение размеров');
|
||
Result:=new Matrix(pa.RowCount,pb.ColCount);
|
||
for var i:=0 to pa.RowCount-1 do
|
||
for var j:=0 to pb.ColCount-1 do begin
|
||
var s:=0.0;
|
||
for var k:=0 to pa.ColCount-1 do
|
||
s+=pa.Value[i,k]*pb.Value[k,j];
|
||
Result.Value[i,j]:=s
|
||
end
|
||
end;
|
||
|
||
class function operator * (pa:Matrix; pb:Vector):Vector;
|
||
// умножение матрицы на вектор
|
||
begin
|
||
Assert(pa.ColCount=pb.Length,
|
||
'Matrix: умножение матрицы на вектор. Несовпадение размеров');
|
||
Result:=new Vector(pa.RowCount);
|
||
for var i:=0 to pa.RowCount-1 do begin
|
||
var s:=0.0;
|
||
for var j:=0 to pa.ColCount-1 do
|
||
s+=pa.Value[i,j]*pb.Value[j];
|
||
Result.Value[i]:=s
|
||
end
|
||
end;
|
||
|
||
class function operator * (pa:Vector; pb:Matrix):Vector;
|
||
// умножение вектора на матрицу
|
||
begin
|
||
Assert(pa.Length=pb.RowCount,
|
||
'Matrix: умножение вектора на матрицу. Несовпадение размеров');
|
||
Result:=new Vector(pb.ColCount);
|
||
for var j:=0 to pb.ColCount-1 do begin
|
||
var s:=0.0;
|
||
for var i:=0 to pb.RowCount-1 do
|
||
s+=pa.Value[i]*pb.Value[i,j];
|
||
Result.Value[j]:=s
|
||
end
|
||
end;
|
||
|
||
/// возвращает диагональную матрицу n*n на базе значения r
|
||
class function Diag(n:integer; r:real):Matrix;
|
||
begin
|
||
Result:=new Matrix(n,n);
|
||
for var i:=0 to n-1 do
|
||
for var j:=0 to n-1 do
|
||
Result.Value[i,j]:=i=j?r:0.0
|
||
end;
|
||
|
||
/// выделение строки
|
||
function Row(k:integer; base:integer:=1):Vector;
|
||
begin
|
||
if base<>0 then base:=1;
|
||
Result:=new Vector(ColCount);
|
||
Result.Value:=Value.Row(k-base)
|
||
end;
|
||
|
||
/// замена строки
|
||
procedure SetRow(v:Vector; k:integer; base:integer:=1);
|
||
begin
|
||
if base<>0 then base:=1;
|
||
Value.SetRow(k-base,v.Value)
|
||
end;
|
||
|
||
/// вставка строки послк указанной
|
||
procedure InsertRowAfter(v:Vector; k:integer; base:integer:=1);
|
||
begin
|
||
if base<>0 then base:=1;
|
||
var k1:=k-base;
|
||
Assert((k1 >= 0) and (k1 < RowCount),
|
||
'Matrix.InsertRowAfter: Недопустимый номер строки');
|
||
Assert(ColCount=v.Length,
|
||
'Matrix.InsertRowAfter: Несовместимая длина вектора');
|
||
var tm:=new Matrix(RowCount+1,ColCount);
|
||
for var i:=0 to k1 do
|
||
for var j:=0 to ColCount-1 do tm.Value[i,j]:=Value[i,j];
|
||
for var j:=0 to ColCount-1 do tm.Value[k1+1,j]:=v.Value[j];
|
||
for var i:=k1+2 to RowCount do
|
||
for var j:=0 to ColCount-1 do tm.Value[i,j]:=Value[i-1,j];
|
||
RowCount+=1;
|
||
Value:=new real[RowCount,ColCount];
|
||
for var i:=0 to RowCount-1 do
|
||
for var j:=0 to ColCount-1 do Value[i,j]:=tm.Value[i,j];
|
||
end;
|
||
|
||
/// вставка строки перед указанной
|
||
procedure InsertRowBefore(v:Vector; k:integer; base:integer:=1);
|
||
begin
|
||
if base<>0 then base:=1;
|
||
var k1:=k-base;
|
||
Assert((k1 >= 0) and (k1 < RowCount),
|
||
'Matrix.InsertRowBefore: Недопустимый номер строки');
|
||
Assert(ColCount=v.Length,
|
||
'Matrix.InsertRowBefore: Несовместимая длина вектора');
|
||
var tm:=new Matrix(RowCount+1,ColCount);
|
||
for var i:=0 to k1-1 do
|
||
for var j:=0 to ColCount-1 do tm.Value[i,j]:=Value[i,j];
|
||
for var j:=0 to ColCount-1 do tm.Value[k1,j]:=v.Value[j];
|
||
for var i:=k1+1 to RowCount do
|
||
for var j:=0 to ColCount-1 do tm.Value[i,j]:=Value[i-1,j];
|
||
RowCount+=1;
|
||
Value:=new real[RowCount,ColCount];
|
||
for var i:=0 to RowCount-1 do
|
||
for var j:=0 to ColCount-1 do Value[i,j]:=tm.Value[i,j];
|
||
end;
|
||
|
||
/// удаление строки
|
||
procedure DeleteRow(k:integer; base:integer:=1);
|
||
begin
|
||
if base<>0 then base:=1;
|
||
var k1:=k-base;
|
||
Assert((k1 >= 0) and (k1 < RowCount),
|
||
'Matrix.DeleteRowAfter: Недопустимый номер строки');
|
||
var tm:=new Matrix(RowCount-1,ColCount);
|
||
for var i:=0 to k1-1 do
|
||
for var j:=0 to ColCount-1 do tm.Value[i,j]:=Value[i,j];
|
||
for var i:=k1+1 to RowCount-1 do
|
||
for var j:=0 to ColCount-1 do tm.Value[i-1,j]:=Value[i,j];
|
||
RowCount-=1;
|
||
Value:=new real[RowCount,ColCount];
|
||
for var i:=0 to RowCount-1 do
|
||
for var j:=0 to ColCount-1 do Value[i,j]:=tm.Value[i,j];
|
||
end;
|
||
|
||
/// Суммирование j-й строки с i-й, результат в i-й строке
|
||
procedure AddRow(i,j:integer; base:integer:=1);
|
||
begin
|
||
if base<>0 then base:=1;
|
||
for var k:=0 to ColCount-1 do
|
||
Value[i-base,k]+=Value[j-base,k]
|
||
end;
|
||
|
||
/// Вычитание j-й строки из i-й, результат в i-й строке
|
||
procedure SubRow(i,j:integer; base:integer:=1);
|
||
begin
|
||
if base<>0 then base:=1;
|
||
for var k:=0 to ColCount-1 do
|
||
Value[i-base,k]-=Value[j-base,k]
|
||
end;
|
||
|
||
/// Умножение строки на скаляр
|
||
procedure MultRow(k:integer; r:real; base:integer:=1);
|
||
begin
|
||
if base<>0 then base:=1;
|
||
for var j:=0 to ColCount-1 do
|
||
Value[k-base,j]*=r
|
||
end;
|
||
|
||
/// Обмен местами строк i и j
|
||
procedure SwapRows(i,j:integer; base:integer:=1);
|
||
begin
|
||
if base<>0 then base:=1;
|
||
Value.SwapRows(i-base,j-base)
|
||
end;
|
||
|
||
/// выделение столбца
|
||
function Col(k:integer; base:integer:=1):Vector;
|
||
begin
|
||
if base<>0 then base:=1;
|
||
Result:=new Vector(RowCount);
|
||
Result.Value:=Value.Col(k-base)
|
||
end;
|
||
|
||
/// замена столбца
|
||
procedure SetCol(v:Vector; k:integer; base:integer:=1);
|
||
begin
|
||
if base<>0 then base:=1;
|
||
Value.SetCol(k-base,v.Value)
|
||
end;
|
||
|
||
/// вставка столбца послк указанного
|
||
procedure InsertColAfter(v:Vector; k:integer; base:integer:=1);
|
||
begin
|
||
if base<>0 then base:=1;
|
||
var k1:=k-base;
|
||
Assert((k1 >= 0) and (k1 < ColCount),
|
||
'Matrix.InsertColAfter: Недопустимый номер столбца');
|
||
Assert(RowCount=v.Length,
|
||
'Matrix.InsertColAfter: Несовместимая длина вектора');
|
||
var tm:=new Matrix(RowCount,ColCount+1);
|
||
for var j:=0 to k1 do
|
||
for var i:=0 to RowCount-1 do tm.Value[i,j]:=Value[i,j];
|
||
for var i:=0 to RowCount-1 do tm.Value[i,k1+1]:=v.Value[i];
|
||
for var j:=k1+2 to ColCount do
|
||
for var i:=0 to RowCount-1 do tm.Value[i,j]:=Value[i,j-1];
|
||
ColCount+=1;
|
||
Value:=new real[RowCount,ColCount];
|
||
for var i:=0 to RowCount-1 do
|
||
for var j:=0 to ColCount-1 do Value[i,j]:=tm.Value[i,j];
|
||
end;
|
||
|
||
/// вставка столбца перед указанным
|
||
procedure InsertColBefore(v:Vector; k:integer; base:integer:=1);
|
||
begin
|
||
if base<>0 then base:=1;
|
||
var k1:=k-base;
|
||
Assert((k1 >= 0) and (k1 < ColCount),
|
||
'Matrix.InsertColBefore: Недопустимый номер столбца');
|
||
Assert(RowCount=v.Length,
|
||
'Matrix.InsertColBefore: Несовместимая длина вектора');
|
||
var tm:=new Matrix(RowCount,ColCount+1);
|
||
for var j:=0 to k1-1 do
|
||
for var i:=0 to RowCount-1 do tm.Value[i,j]:=Value[i,j];
|
||
for var i:=0 to RowCount-1 do tm.Value[i,k1]:=v.Value[i];
|
||
for var j:=k1+1 to ColCount do
|
||
for var i:=0 to RowCount-1 do tm.Value[i,j]:=Value[i,j-1];
|
||
ColCount+=1;
|
||
Value:=new real[RowCount,ColCount];
|
||
for var i:=0 to RowCount-1 do
|
||
for var j:=0 to ColCount-1 do Value[i,j]:=tm.Value[i,j];
|
||
end;
|
||
|
||
/// удаление столбца
|
||
procedure DeleteCol(k:integer; base:integer:=1);
|
||
begin
|
||
if base<>0 then base:=1;
|
||
var k1:=k-base;
|
||
Assert((k1 >= 0) and (k1 < ColCount),
|
||
'Matrix.DeleteCol: Недопустимый номер столбца');
|
||
var tm:=new Matrix(RowCount,ColCount-1);
|
||
for var j:=0 to k1-1 do
|
||
for var i:=0 to RowCount-1 do tm.Value[i,j]:=Value[i,j];
|
||
for var j:=k1+1 to ColCount-1 do
|
||
for var i:=0 to RowCount-1 do tm.Value[i,j-1]:=Value[i,j];
|
||
ColCount-=1;
|
||
Value:=new real[RowCount,ColCount];
|
||
for var i:=0 to RowCount-1 do
|
||
for var j:=0 to ColCount-1 do Value[i,j]:=tm.Value[i,j];
|
||
end;
|
||
|
||
/// Суммирование j-го столбца с i-м, результат в i-м столбце
|
||
procedure AddCol(i,j:integer; base:integer:=1);
|
||
begin
|
||
if base<>0 then base:=1;
|
||
for var k:=0 to RowCount-1 do
|
||
Value[k,i-base]+=Value[k,j-base]
|
||
end;
|
||
|
||
/// Вычитание j-го столбца из i-го, результат в i-м столбце
|
||
procedure SubCol(i,j:integer; base:integer:=1);
|
||
begin
|
||
if base<>0 then base:=1;
|
||
for var k:=0 to RowCount-1 do
|
||
Value[k,i-base]-=Value[k,j-base]
|
||
end;
|
||
|
||
/// Умножение столбца на скаляр
|
||
procedure MultCol(k:integer; r:real; base:integer:=1);
|
||
begin
|
||
if base<>0 then base:=1;
|
||
for var i:=0 to RowCount-1 do
|
||
Value[i,k-base]*=r
|
||
end;
|
||
|
||
/// Обмен местами столбцов i и j
|
||
procedure SwapCols(i,j:integer; base:integer:=1);
|
||
begin
|
||
if base<>0 then base:=1;
|
||
Value.SwapCols(i-base,j-base)
|
||
end;
|
||
|
||
/// Транспонирование матрицы
|
||
function Transpose:Matrix;
|
||
begin
|
||
Result:=new Matrix(ColCount,RowCount);
|
||
Result.Value:=PABCSystem.Transpose(Value);
|
||
end;
|
||
|
||
/// вывод матрицы на экран
|
||
function Println(w:integer:=6; d:integer:=2):Matrix;
|
||
begin
|
||
Result:=new Matrix(RowCount,ColCount);
|
||
for var i:=0 to RowCount-1 do begin
|
||
for var j:=0 to ColCount-1 do begin
|
||
Result.Value[i,j]:=Value[i,j];
|
||
Write(Result.Value[i,j]:w:d)
|
||
end;
|
||
Writeln
|
||
end
|
||
end;
|
||
|
||
/// детерминант матрицы
|
||
function Det:real;
|
||
begin
|
||
var r:=new Decomp(Value);
|
||
Result:=r.det
|
||
end;
|
||
|
||
/// обращение квадратной матрицы
|
||
function Inv:Matrix;
|
||
begin
|
||
Assert(RowCount=ColCount,'Matrix.Inv: матрица не квадратная');
|
||
Result:=new Matrix(RowCount,ColCount);
|
||
var r:=new Decomp(Value);
|
||
if r.cond<MaxReal then begin
|
||
for var j:=0 to ColCount-1 do begin
|
||
var v:=ArrGen(RowCount,i->i=j?1.0:0.0);
|
||
r.Solve(v);
|
||
for var i:=0 to RowCount-1 do Result.Value[i,j]:=v[i]
|
||
end
|
||
end
|
||
else Result.Value:=MatrFill(RowCount,ColCount,MaxReal);
|
||
end;
|
||
|
||
/// Решение системы линейных алгебраических уравнения (СЛАУ)
|
||
function SLAU(vb:Vector; var cond:real):Vector;
|
||
// Возвращает вектор, содержащий решение СЛАУ.
|
||
// Заносит в cond число обусловленности исходной матрицы.
|
||
// Если cond велико, исходная матрица вырождена или близка к этому и
|
||
// вектором результата пользоваться нельзя.
|
||
// Является оболочкой к самостоятельному классу Decomp.
|
||
begin
|
||
Assert(vb.Length=ColCount,
|
||
'Matrix.SLAU: длина вектора vb несовместима с размером матрицы');
|
||
var r:=new Decomp(Value);
|
||
cond:=r.cond;
|
||
Result:=new Vector(vb.Length);
|
||
Result.Value:=PABCSystem.Copy(vb.Value);
|
||
r.Solve(Result.Value);
|
||
end;
|
||
|
||
/// возвращает копию матрицы
|
||
function Copy:Matrix;
|
||
begin
|
||
Result:=new Matrix(RowCount,ColCount);
|
||
for var i:=0 to RowCount-1 do
|
||
for var j:=0 to ColCount-1 do
|
||
Result.Value[i,j]:=Value[i,j]
|
||
end;
|
||
|
||
end;
|
||
{$endregion}
|
||
|
||
{$region Polynom}
|
||
type
|
||
/// Полиномиальная арифметика
|
||
Polynom=class
|
||
|
||
public
|
||
a:array of real; // коэффициенты по возрастанию степени
|
||
n:integer; // кол-во элементов в полиноме
|
||
eps:real; // оценка точности при экономизации полинома
|
||
|
||
constructor(params p:array of real);
|
||
begin
|
||
a:=PABCSystem.Copy(p);
|
||
n:=a.Length
|
||
end;
|
||
|
||
constructor Create(n:integer);
|
||
begin
|
||
a:=ArrFill(n,0.0);
|
||
Self.n:=n
|
||
end;
|
||
|
||
/// Значение полинома P(x)
|
||
function Value(x:real):real;
|
||
begin
|
||
Result:=a[0];
|
||
var k:=1.0;
|
||
for var i:=1 to n-1 do begin
|
||
k*=x;
|
||
Result+=a[i]*k
|
||
end;
|
||
end;
|
||
|
||
/// коэффициенты производной полинома
|
||
function PDif:Polynom;
|
||
begin
|
||
Result:=new Polynom(n-1);
|
||
for var i:=0 to n-2 do Result.a[i]:=a[i+1]*(i+1)
|
||
end;
|
||
|
||
/// коэффициенты неопределенного интеграла от полинома
|
||
function PInt:Polynom;
|
||
begin
|
||
Result:=new Polynom(Arr(0.0)+a);
|
||
for var i:=1 to Result.n-1 do Result.a[i]/=i
|
||
end;
|
||
|
||
class function operator + (p,q:Polynom):Polynom;
|
||
begin
|
||
var m:=Max(p.n,q.n);
|
||
Result:=new Polynom(m);
|
||
if p.n>q.n then begin
|
||
for var i:=0 to q.n-1 do Result.a[i]:=p.a[i]+q.a[i];
|
||
for var i:=q.n to m-1 do Result.a[i]:=p.a[i]
|
||
end
|
||
else begin
|
||
for var i:=0 to p.n-1 do Result.a[i]:=p.a[i]+q.a[i];
|
||
for var i:=p.n to m-1 do Result.a[i]:=q.a[i]
|
||
end
|
||
end;
|
||
|
||
class function operator + (r:real; p:Polynom):Polynom;
|
||
begin
|
||
Result:=new Polynom(p.a);
|
||
Result.a[0]+=r
|
||
end;
|
||
|
||
class function operator + (p:Polynom; r:real):=r+p;
|
||
|
||
class function operator - (p:Polynom):=
|
||
new Polynom(p.a.Select(x->-x).ToArray);
|
||
|
||
class function operator - (p,q:Polynom):Polynom;
|
||
begin
|
||
var m:=Max(p.n,q.n);
|
||
Result:=new Polynom(m);
|
||
if p.n>q.n then begin
|
||
for var i:=0 to q.n-1 do Result.a[i]:=p.a[i]-q.a[i];
|
||
for var i:=q.n to m-1 do Result.a[i]:=p.a[i]
|
||
end
|
||
else begin
|
||
for var i:=0 to p.n-1 do Result.a[i]:=p.a[i]-q.a[i];
|
||
for var i:=p.n to m-1 do Result.a[i]:=-q.a[i]
|
||
end
|
||
end;
|
||
|
||
class function operator - (r:real; p:Polynom):Polynom;
|
||
begin
|
||
Result:=-p;
|
||
Result.a[0]+=r
|
||
end;
|
||
|
||
class function operator - (p:Polynom; r:real):Polynom;
|
||
begin
|
||
Result:=new Polynom(p.a);
|
||
Result.a[0]-=r
|
||
end;
|
||
|
||
class function operator * (p,q:Polynom):Polynom;
|
||
begin
|
||
var m:=p.n+q.n-1;
|
||
Result:=new Polynom(m);
|
||
for var i:=0 to p.n-1 do
|
||
for var j:=0 to q.n-1 do
|
||
Result.a[i+j]+=p.a[i]*q.a[j]
|
||
end;
|
||
|
||
class function operator * (r:real; p:Polynom):Polynom;
|
||
begin
|
||
Result:=new Polynom(p.a);
|
||
for var i:=0 to p.n-1 do Result.a[i]*=r
|
||
end;
|
||
|
||
class function operator * (p:Polynom; r:real):=r*p;
|
||
|
||
class function operator = (p,q:Polynom):boolean;
|
||
begin
|
||
if p.n<>q.n then Result:=false
|
||
else begin
|
||
Result:=true;
|
||
for var i:=0 to p.n-1 do
|
||
if Abs(p.a[i]-q.a[i])>0 then begin
|
||
Result:=false; Exit
|
||
end
|
||
end
|
||
end;
|
||
|
||
/// Копия полинома
|
||
function Copy:=new Polynom(a);
|
||
|
||
class function operator / (p,q:Polynom):(Polynom,Polynom);
|
||
begin
|
||
var np:=p.n-q.n+1;
|
||
if np<=0 then begin
|
||
Result:=(new Polynom(1),p.Copy);
|
||
Exit
|
||
end;
|
||
var pp:=new Polynom(np);
|
||
var s:=p.Copy;
|
||
var ns:=q.n-1;
|
||
var i:=np;
|
||
repeat
|
||
pp.a[i-1]:=s.a[i+ns-1]/q.a.Last;
|
||
for var k:=0 to ns-1 do begin
|
||
var j:=k-1+i;
|
||
s.a[j]-=pp.a[i-1]*q.a[k]
|
||
end;
|
||
i-=1
|
||
until i<=0;
|
||
SetLength(s.a,ns);
|
||
s.n:=ns;
|
||
Result:=(pp,s)
|
||
end;
|
||
|
||
class function operator / (p:Polynom; q:real):Polynom;
|
||
begin
|
||
Result:=p.Copy;
|
||
for var i:=0 to Result.n-1 do Result.a[i]/=q
|
||
end;
|
||
|
||
class function operator / (p:real; q:Polynom):=(p,q.Copy);
|
||
|
||
/// Вывод коэффициентов полинома с разделителем
|
||
function Print(delim:char:=' '):Polynom;
|
||
begin
|
||
a.Print(delim);
|
||
Result:=Self.Copy
|
||
end;
|
||
|
||
/// Замена аргумента х на at+b
|
||
function Polyx(a,b:real):Polynom;
|
||
begin
|
||
Result:=new Polynom(n);
|
||
var z:=new real[n];
|
||
var w:=ArrFill(n,1.0);
|
||
Result.a[0]:=Self.a[0];
|
||
z[0]:=1;
|
||
for var i:=1 to n-1 do begin
|
||
z[i]:=b*z[i-1];
|
||
Result.a[0]+=Self.a[i]*z[i]
|
||
end;
|
||
for var j:=1 to n-1 do begin
|
||
w[0]*=a; Result.a[j]:=Self.a[j]*w[0];
|
||
for var i:=j+1 to n-1 do begin
|
||
var k:=i-j;
|
||
w[k]:=a*w[k]+w[k-1];
|
||
Result.a[j]+=Self.a[i]*w[k]*z[k]
|
||
end
|
||
end
|
||
end;
|
||
|
||
/// Экономизация полинома на [-1;1]
|
||
function EconomSym(h,limit:real):Polynom;
|
||
// Переработка алгоритма 38б
|
||
// В кн.: Агеев М.И., Алик В.П., Марков Ю.И.
|
||
// Библиотека алгоритмов 1б-50б. (Справочное пособие).
|
||
// М., "Сов.радио", 1975
|
||
begin
|
||
var n1:=n-1;
|
||
var t:=new real[n];
|
||
Result:=Self.Copy;
|
||
var s:real;
|
||
eps:=0;
|
||
while n1>=2 do begin
|
||
t[n1]:=-Result.a[n1];
|
||
var k:=n1;
|
||
while k>=2 do begin
|
||
t[k-2]:=-t[k]*Sqr(h)*k*(k-1)/((n1+k-2)*(n1-k+2));
|
||
k-=2
|
||
end;
|
||
s:=n1.IsEven?a[0]:t[1]/n1;
|
||
if eps+Abs(s)<limit then begin
|
||
eps+=Abs(s);
|
||
k:=n1;
|
||
while k>=0 do begin Result.a[k]+=t[k]; k-=2 end;
|
||
n1-=1;
|
||
end
|
||
else break
|
||
end;
|
||
SetLength(Result.a,n1+1);
|
||
Result.n:=n1+1;
|
||
Result.eps:=eps
|
||
end;
|
||
|
||
/// Экономизация полинома на [0,1]
|
||
function EconomUnsym(h,limit:real):Polynom;
|
||
// Переработка алгоритма 37б
|
||
// В кн.: Агеев М.И., Алик В.П., Марков Ю.И.
|
||
// Библиотека алгоритмов 1б-50б. (Справочное пособие).
|
||
// М., "Сов.радио", 1975
|
||
begin
|
||
var n1:=n-1;
|
||
var t:=new real[n];
|
||
Result:=Self.Copy;
|
||
eps:=0.0;
|
||
while n1>=1 do begin
|
||
t[n1]:=-Result.a[n1];
|
||
for var k:=n1 downto 1 do
|
||
t[k-1]:=-t[k]*h*k*(k-0.5)/((n1+k-1)*(n1-k+1));
|
||
if eps+Abs(t[0])<limit then begin
|
||
eps+=Abs(t[0]);
|
||
for var k:=n1 downto 0 do Result.a[k]+=t[k];
|
||
n1-=1;
|
||
end
|
||
else break
|
||
end;
|
||
SetLength(Result.a,n1+1);
|
||
Result.n:=n1+1;
|
||
Result.eps:=eps
|
||
end;
|
||
|
||
/// Вывод коэффициентов полинома с разделителем и сменой строки
|
||
function Println(delim:char:=' '):Polynom;
|
||
begin
|
||
a.Println(delim);
|
||
Result:=Self.Copy
|
||
end;
|
||
|
||
/// Вывод полинома в "математическом" виде со сменой строки
|
||
procedure PrintlnBeauty(p:char:='x');
|
||
begin
|
||
var s:string:='';
|
||
for var i:=Self.n-1 downto 0 do begin
|
||
var t:=Self.a[i];
|
||
case i of
|
||
0: if t<0 then s+=t.ToString
|
||
else if t>0 then s+='+'+t.ToString
|
||
else if s.Length=0 then s:='0';
|
||
1: if t<0 then
|
||
if t<>-1 then s+=t.ToString+p
|
||
else s+='-'+p
|
||
else if t>0 then
|
||
if t<>1 then s+='+'+t.ToString+p
|
||
else s+='+'+p;
|
||
else
|
||
if t<0 then
|
||
if t<>-1 then s+=t.ToString+p+'^'+i.ToString
|
||
else s+='-'+p+'^'+i.ToString
|
||
else if t>0 then
|
||
if t<>1 then s+='+'+t.ToString+p+'^'+i.ToString
|
||
else s+='+'+p+'^'+i.ToString
|
||
end
|
||
end;
|
||
if s[1]='+' then s:=s[2:];
|
||
Writeln(s)
|
||
end;
|
||
|
||
end;
|
||
{$endregion}
|
||
|
||
{$region PolRt}
|
||
type
|
||
/// Нахождение всех корней полинома с действительными коэффициентами
|
||
/// методом Ньютона-Рафсона (не выше 36-й степени)
|
||
PolRt=class
|
||
|
||
private
|
||
p:Polynom;
|
||
|
||
public
|
||
ier:integer;
|
||
|
||
constructor(p:Polynom);
|
||
begin
|
||
Self.p:=p;
|
||
end;
|
||
|
||
/// Вывод корней полинома в массив
|
||
function Value:array of complex;
|
||
//
|
||
// Существенно переработанная для языка PascalABC.NET 3.2
|
||
// подпрограмма DPOLRT из библиотеки SSP
|
||
// (она же - "Пакет научных программ на языке Фортран-4").
|
||
//
|
||
// Коэффициенты полинома задаются в порядке возрастания степени.
|
||
// ier=0 - ошибок не обнаружено
|
||
// ier=1 - мало элементов для построения полинома
|
||
// ier=2 - степень полинома превышает 36
|
||
// ier=3 - не удалось достигнуть приемлемой точности решения за 500 шагов
|
||
// ier=4 - коэффициент при старшей степени полинома нулевой
|
||
//
|
||
begin
|
||
var m:=p.n-2;
|
||
if m>35 then begin ier:=2; Exit end;
|
||
if p.a[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:=p.a.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(p.a[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;
|
||
|
||
end;
|
||
{$endregion}
|
||
|
||
{$region Quanc8}
|
||
type
|
||
/// Адаптивная квадратурная программа численного интегрирования
|
||
Quanc8=class
|
||
|
||
private
|
||
|
||
fun:real->real;
|
||
a,b,abserr,relerr:real;
|
||
|
||
public
|
||
|
||
constructor(f:real->real; a,b,abserr,relerr:real);
|
||
begin
|
||
fun:=f;
|
||
Self.a:=a;
|
||
Self.b:=b;
|
||
Self.abserr:=abserr;
|
||
Self.relerr:=relerr;
|
||
end;
|
||
|
||
/// Вычисление и возврат значения интеграла в виде кортежа
|
||
/// Непосредственно результат первый в кортеже
|
||
function Value:(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;
|
||
|
||
end;
|
||
{$endregion}
|
||
|
||
{$region RKF45}
|
||
type
|
||
/// Метод Рунге-Кутты-Фельберга 4-5 порядка
|
||
RKF45=class
|
||
|
||
// За основу взят текст программы RKF45 на языке Fortran-90,
|
||
// находящейся на сайте https://www.sc.fsu.edu
|
||
// ( Dept. of Scientific Computing, Florida State University ),
|
||
// распространяемой на условиях GNU LGPL
|
||
|
||
private
|
||
|
||
neqn:integer; // количество уравнений системы
|
||
abserr_save:real;
|
||
relerr_save:real;
|
||
h:real;
|
||
epsilon:=1.2e-15;
|
||
flag_save:integer;
|
||
init:integer;
|
||
kflag:integer;
|
||
kop:integer;
|
||
nfe:integer;
|
||
f1:array of real;
|
||
f2:array of real;
|
||
f3:array of real;
|
||
f4:array of real;
|
||
f5:array of real;
|
||
p:procedure(x:real; y,yp:array of real);
|
||
|
||
function Sign(a,b:real):=b<0?-a:a;
|
||
|
||
procedure Fehl(s:array of real; t,h:real);
|
||
//
|
||
// Метод Рунге-Кутты-Фельберга 4-5 порядка.
|
||
// Интегрируется система обыкновенных дифференциальных уравнений 1-го порядка
|
||
// вида dy(i)/dt = f(t, y(1), y(2) ... ,
|
||
// где значения функции y(i) и значения её производных yp(i) заданы
|
||
// в начальной точке t.
|
||
// Метод продолжает решение на фиксированный шаг h и помещает в массив s
|
||
// приближение к решению в точке t+h, имеющее пятый порядок точности.
|
||
// В формулах произведена группировка с целью уменьшить потерю
|
||
// верных знаков.
|
||
// Чтобы можно было различать значения независимых аргументов, при вызове
|
||
// Fehl не следует задавать для h значение, меньшее ошибки округления t,
|
||
// умноженной на 13.
|
||
//
|
||
begin
|
||
var ch:=h/4.0;
|
||
for var i:=0 to neqn-1 do f5[i]:=y[i]+ch*yp[i];
|
||
p(t+ch,f5,f1); ch:=3.0*h/32.0;
|
||
for var i:=0 to neqn-1 do f5[i]:=y[i]+ch*(yp[i]+3.0*f1[i]);
|
||
p(t+3.0*h/8.0,f5,f2); ch:=h/2197.0;
|
||
for var i:=0 to neqn-1 do
|
||
f5[i]:=y[i]+ch*(1932.0*yp[i]+7296.0*f2[i]-7200.0*f1[i]);
|
||
p(t+12.0*h/13.0,f5,f3); ch:=h/4104.0;
|
||
for var i:=0 to neqn-1 do
|
||
f5[i]:=y[i]+ch*(8341.0*yp[i]-845.0*f3[i]+29440.0*f2[i]-32832.0*f1[i]);
|
||
p(t+h,f5,f4); ch:=h/20520.0;
|
||
for var i:=0 to neqn-1 do
|
||
f1[i]:=y[i]+ch*(-6080.0*yp[i]+9295.0*f3[i]-5643.0*f4[i]+41040.0*f1[i]-
|
||
28352.0*f2[i]);
|
||
p(t+h/2.0,f1,f5); ch:=h/7618050.0;
|
||
for var i:=0 to neqn-1 do
|
||
s[i]:=y[i]+ch*(902880.0*yp[i]+3855735.0*f3[i]-1371249.0*f4[i]+
|
||
3953664.0*f2[i]+277020.0*f5[i])
|
||
end;
|
||
|
||
procedure FlagDev;
|
||
// обработчик состояния флага завершенеия шага flag
|
||
// Выводит сообщения при MsgOn=true
|
||
begin
|
||
if MsgOn then
|
||
case flag of
|
||
-5:Writeln('RKF45: Критическая ошибка',NewLine,'kFlag=5 и AbsErr=0');
|
||
-4:Writeln('RKF45: Критическая ошибка',NewLine,
|
||
'kFlag=6 и RelErr<=RelErr_save и AbsErr<=AbsErr_save');
|
||
-3:Writeln('RKF45: Критическая ошибка',NewLine,
|
||
'Невозможно продолжить интегрирование',NewLine,
|
||
'Пользователь не отреагировал на flag=5..8');
|
||
1: Writeln('Ненадлежащий вызов');
|
||
3: Writeln('Установлено: relerr,abserr = ',relerr,',',abserr);
|
||
4: Writeln('Требуется слишком много шагов');
|
||
5: begin
|
||
abserr:=0.2e-5;
|
||
Writeln('Установлено: relerr,abserr = ',relerr,',',abserr);
|
||
flag:=0
|
||
end;
|
||
6: begin
|
||
relerr*=10.0;
|
||
Writeln('Установлено: relerr,abserr = ',relerr,',',abserr);
|
||
flag:=0
|
||
end;
|
||
7: begin
|
||
Writeln('Сложно получить требуемую точность');
|
||
flag:=0
|
||
end;
|
||
8: Writeln('Ненадлежащий вызов')
|
||
end;
|
||
end;
|
||
|
||
public
|
||
|
||
flag:integer;
|
||
relerr:real;
|
||
abserr:real;
|
||
y:array of real; // функция
|
||
yp:array of real; // первая производная
|
||
MsgOn:boolean;
|
||
|
||
constructor (f:procedure(x:real; y,yp:array of real);
|
||
y:array of real; abserr,relerr:real; MsgOn:boolean:=true);
|
||
begin
|
||
neqn:=y.Length;
|
||
Self.y:=Copy(y);
|
||
yp:=new real[neqn];
|
||
f1:=new real[neqn];
|
||
f2:=new real[neqn];
|
||
f3:=new real[neqn];
|
||
f4:=new real[neqn];
|
||
f5:=new real[neqn];
|
||
Self.epsilon:=epsilon;
|
||
flag:=1;
|
||
Self.abserr:=abserr;
|
||
Self.relerr:=relerr;
|
||
Self.MsgOn:=MsgOn;
|
||
(abserr_save,relerr_save,h):=(-1.0,-1.0,-1.0);
|
||
(flag_save,init,kflag,kop,nfe):=(-1000,-1000,-1000,-1,-1);
|
||
p:=f;
|
||
end;
|
||
|
||
procedure Solve(var t:real; tout:real);
|
||
// Реализует метод Runge-Kutta-Fehlberg
|
||
// Обычно данный класс выполняет интегрирование на отрезке
|
||
// от t до t_out, но если требуется одношаговвй интегратор,
|
||
// можно с его помощью сделать небольшой шажок в направлении t_out.
|
||
// После каждого шага значения t м t_out устанавливаются так,
|
||
// чтобы для следующего шага достаточно только еще раз
|
||
// обратиться к методу Solve.
|
||
//
|
||
// Перед первым обращением пользователь должен
|
||
// - обеспечить наличие процедуры с интерфейсом вида F(t,y,yp),
|
||
// вычисляющей при заданном аргументе t значения функции y и её первой
|
||
// производной yp на правом конце интервала;
|
||
// - задать начальные значения параметоров y, t, tout, relerr, abserr и flag.
|
||
// Обычно t задает начальную точку для решения,
|
||
// y определяет начальные условия, а flag полагают равным +1.
|
||
//
|
||
// Cтандартно значение flag перед первым вызовом устанавливыается в 1,
|
||
// а дальнейшее управление берет на себя программа. Flag=-1 задается,
|
||
// если нужен одношаговый режим.
|
||
// После завершения метод Solve возвращает значение Flag=2 в обычном режиме
|
||
// или -2 в одношаговом, что позволяет нормально продолжить работу.
|
||
// (Значения +1 и -1 указывают на необходимость некоторой внутренней
|
||
// инициализации, а 2 и -2 вызывают её обход)
|
||
//
|
||
// Если возвращенный результат равен заданному tout, пользователь
|
||
// должен определить новое значение для следующего шага и продолжить работу.
|
||
// В одношаговом режиме пользователь, получив Flag=-2, должен понимать,
|
||
// что решение - это лишь очередная итерация в направлении tout и следует
|
||
// повторить вызов метода Solve.
|
||
// Значение Flag=2 независимо от режима указывает, что tout достигнуто.
|
||
// Если нужно продолжить работу в одношаговом режиме, следует задать Flag=-2.
|
||
//
|
||
// В некоторых случаях могут возникать ошибки. Понять их причину и
|
||
// принять соответствующее решение помогает значение Flag.
|
||
//
|
||
// 3: интегрирование не завершено, поскольку заданное значение границы
|
||
// льносительной ошибки relerr слишкои мало. Для продолжения необходимо
|
||
// увеличить значение relerr. Если результат устраивает установите
|
||
// Flag=2 и продолжайтк работу.
|
||
//
|
||
// 4: интегрирование не завершено, поскольку потребовалось более 3000
|
||
// вычислений значения производной, что соответствует примерно 500 шагам.
|
||
// Пользователь может повторить обращение, при этом счетчик
|
||
// количества вычислений функции будет сброшен в ноль и, возможно, это
|
||
// даст результат.
|
||
// 5: интегрирование не завершено, поскольку решение обратилось в ноль,
|
||
// вследствие чего тест одной только относительной ошибки не проходит.
|
||
// Пользователь должен задать ненулевое значение abserr и продолжить работу.
|
||
// Возможно, использование одношагового интегратора будет хорошим выходом.
|
||
// 6: интегрирование не завершено, поскольку требуемая точность не может
|
||
// быть достигнута даже при наименьшем допустимом шаге. Пользователь должен
|
||
// увеличить границу допустимой погрешности abserr и/или reller прежде чем
|
||
// пытаться продолжить решение. Также, потребуется установить flag=2 (или -2
|
||
// для вызова одношагового интегратора). flag=6 возникает в случае
|
||
// обнаружения проблемной точки, в которой или решение меняется очень резко,
|
||
// или присутствует сингулярность.
|
||
// 7: по всей видимости RKF45 неэффективна для для данной задачию. Слишком
|
||
// большое значение требуемых для вычисления точек препятствует выбору шага.
|
||
// Можно попытаться использовать режим одношагового интегратора.
|
||
// Если пользователь попытается продолжить решение, задав flag=2, работа
|
||
// будет завершена.
|
||
// 8: Неверные входные параметры. Возникает, если:
|
||
// - нулевое или отрицательне количество уравнений;
|
||
// - t=tout и при этом flag не равен 1;
|
||
// - relerr или abserr отрицательны;
|
||
// - flag имеет значение 0 или не лежит в интервале [-4;8]
|
||
begin
|
||
var maxnfe:=3000;
|
||
var remin:=1.0E-12;
|
||
// remin - это минимальное доа\пустимое знвчение relerr.
|
||
// Попытки получить с помощью RKF45 более высокую точность обычно
|
||
// безуспешны и обходятся очень дорого.
|
||
var ae,dt,ee,eeoet,eps,esttol,et,hmin,relerr_min,s,scale,tol,toln,ypk:real;
|
||
var hfaild,output:boolean;
|
||
eps:=Self.epsilon; // машинная точность
|
||
// контроль входных значений
|
||
if (neqn<1) then begin
|
||
flag:=8; FlagDev; Exit
|
||
end;
|
||
if (relerr<0.0) or (abserr<0.0) then begin
|
||
flag:=8; FlagDev; Exit
|
||
end;
|
||
if (flag=0) or (flag<-2) or (flag>8) then begin
|
||
flag:=8; FlagDev; Exit
|
||
end;
|
||
var mflag:=Abs(flag);
|
||
if mflag<>1 then begin // обход части кода при первом обращении
|
||
if (t=tout) and (kflag<>3) then begin
|
||
flag:=8; FlagDev; Exit
|
||
end;
|
||
if mflag=2 then begin
|
||
if kflag=3 then begin
|
||
flag:=flag_save; mflag:=Abs(flag)
|
||
end
|
||
else if init=0 then flag:=flag_save
|
||
else if kflag=4 then nfe:=0
|
||
else if (kflag=5) and (abserr=0.0) then begin
|
||
flag:=-5; FlagDev; Exit
|
||
end
|
||
else if (kflag=6) and (relerr<=relerr_save) and
|
||
(abserr<=abserr_save) then begin
|
||
flag:=-4; FlagDev; Exit
|
||
end
|
||
end
|
||
else
|
||
if flag=3 then begin
|
||
flag:=flag_save;
|
||
if kflag=3 then mflag:=Abs(flag);
|
||
end
|
||
else if flag=4 then begin
|
||
nfe:=0; flag:=flag_save;
|
||
if kflag=3 then mflag:=Abs(flag)
|
||
end
|
||
else if (flag=5) and (abserr>0.0) then begin
|
||
flag:=flag_save;
|
||
if kflag=3 then mflag:=Abs(flag)
|
||
end
|
||
else begin
|
||
flag:=-3; FlagDev; Exit
|
||
end
|
||
end;
|
||
flag_save:=flag; kflag:=0;
|
||
relerr_save:=relerr; abserr_save:=abserr; // для будущих проверок
|
||
relerr_min:=2.0*epsilon+remin; // ограниечение по точности
|
||
if relerr<relerr_min then begin
|
||
relerr:=relerr_min; flag:=3; kflag:=3;
|
||
FlagDev; Exit
|
||
end;
|
||
dt:=tout-t;
|
||
if mflag=1 then begin
|
||
init:=0; kop:=0;
|
||
p(t,y,yp);
|
||
nfe:=1;
|
||
if t=tout then begin
|
||
flag:=2; FlagDev; Exit
|
||
end
|
||
end;
|
||
if init=0 then begin
|
||
init:=1; h:=Abs(dt); toln:=0.0;
|
||
for var k:=0 to neqn-1 do begin
|
||
tol:=relerr*Abs(y[k])+abserr;
|
||
if tol>0.0 then begin
|
||
toln:=tol; ypk:=Abs(yp[k]);
|
||
if tol<ypk*Power(h,5) then h:=Power(tol/ypk,0.2)
|
||
end
|
||
end;
|
||
if toln<=0.0 then h:=0.0;
|
||
h:=Max(h,26.0*eps*Max(Abs(t),Abs(dt)));
|
||
flag_save:=flag<0?-2:2
|
||
end;
|
||
if dt<0 then h:=-h;
|
||
if 2.0*Abs(dt)<=Abs(h) then kop+=1;
|
||
if kop=10000 then begin kop:=0; flag:=7; FlagDev; Exit end;
|
||
if Abs(dt)<=26.0*eps*Abs(t) then begin
|
||
t:=tout;
|
||
for var i:=0 to neqn-1 do y[i]+=dt*yp[i];
|
||
p(t,y,yp);
|
||
nfe+=1; flag:=2;
|
||
FlagDev; Exit
|
||
end;
|
||
// интегрирование шаг за шагом
|
||
output:=false; scale:=2.0/relerr; ae:=scale*abserr;
|
||
repeat
|
||
hfaild:=false; hmin:=26.0*eps*Abs(t); dt:=tout-t;
|
||
if 2.0*Abs(h)>Abs(dt) then begin
|
||
if Abs(dt)<=Abs(h) then begin output:=true; h:=dt end
|
||
else h:=0.5*dt
|
||
end;
|
||
// внутренний одношаговый интегратор
|
||
while true do begin
|
||
if maxnfe<nfe then begin (flag,kflag):=(4,4); FlagDev; Exit end;
|
||
Fehl(f1,t,h);
|
||
nfe+=5; eeoet:=0.0;
|
||
for var k:=0 to neqn-1 do begin
|
||
et:=Abs(y[k])+Abs(f1[k])+ae;
|
||
if et<=0.0 then begin flag:=5; FlagDev; Exit end;
|
||
ee:=Abs((-2090.0*yp[k]+(21970.0*f3[k]-15048.0*f4[k]))+
|
||
(22528.0*f2[k]-27360.0*f5[k]));
|
||
eeoet:=Max(eeoet,ee/et)
|
||
end;
|
||
esttol:=Abs(h)*eeoet*scale/752400.0;
|
||
if esttol<=1.0 then Break;
|
||
hfaild:=true; output:=false;
|
||
if esttol<59049.0 then s:=0.9/Power(esttol,0.2)
|
||
else s:=0.1;
|
||
h*=s;
|
||
if Abs(h)<hmin then begin (flag,kflag):=(6,6); FlagDev; Exit end;
|
||
end;
|
||
t+=h; y:=Copy(f1);
|
||
p(t,y,yp);
|
||
nfe+=1;
|
||
if esttol>0.0001889568 then s:=0.9/Power(esttol,0.2)
|
||
else s:=5.0;
|
||
if hfaild then s:=Min(s,1.0);
|
||
h:=Sign(Max(s*Abs(h),hmin),h);
|
||
// окончание кода внутреннего одношагового интегратора
|
||
if output then begin t:=tout; flag:=2; FlagDev; Exit end;
|
||
until flag<=0;
|
||
flag:=-2;
|
||
FlagDev
|
||
end;
|
||
|
||
end;
|
||
{$endregion}
|
||
|
||
{$region RootsIsolation}
|
||
type
|
||
/// Изоляция корней уравнения y(x)=0 на заданном интервале
|
||
RootsIsolation=class
|
||
|
||
private
|
||
f:real->real;
|
||
a,b,h:real;
|
||
|
||
public
|
||
|
||
constructor(f:real->real; a,b,h:real);
|
||
begin
|
||
Self.f:=f;
|
||
Self.a:=a;
|
||
Self.b:=b;
|
||
Self.h:=h
|
||
end;
|
||
|
||
/// Полученные результаты
|
||
function Value:array of (real,real);
|
||
//
|
||
// f - функция, для которой ищутся интервалы изоляции;
|
||
// a,b - границы интервала просмотра [a;b]
|
||
// h - шаг просмотра
|
||
// Функция возвращает массив кортежей (x1,x2), таких что f(x1) и f(x2) имеют
|
||
// разные знаки.
|
||
//
|
||
begin
|
||
var r:=new List<(real,real)>;
|
||
var h2:=h/2;
|
||
var eps:=1e-7;
|
||
var (xa,ya,bk):=(a,f(a),b+h2);
|
||
if Abs(ya)<eps then begin
|
||
r.Add((xa-h2,xa+h2));
|
||
xa+=h2; ya:=f(xa)
|
||
end;
|
||
var xb:=xa+h;
|
||
while xb<=bk do begin
|
||
var yb:=f(xb);
|
||
if Abs(yb)<eps then begin
|
||
r.Add((xb-h2,xb+h2));
|
||
xa:=xb+h2; ya:=f(xa);
|
||
xb:=xa+h
|
||
end
|
||
else if ya*yb<0 then begin
|
||
r.Add((xa,xb));
|
||
(xa,ya):=(xb,yb);
|
||
xb:=xa+h
|
||
end
|
||
else begin
|
||
(xa,ya):=(xb,yb);
|
||
xb:=xa+h
|
||
end
|
||
end;
|
||
Result:=r.ToArray
|
||
end;
|
||
|
||
end;
|
||
{$endregion}
|
||
|
||
{$region Spline}
|
||
type
|
||
/// Интерполяция табличной функции кубическим сплайном
|
||
Spline=class
|
||
|
||
private
|
||
|
||
n:integer; // количество узлов интерполяции
|
||
aK:array[,] of real; // коэффициенты сплайна
|
||
P:array of Point; // узлы интерполяции
|
||
|
||
///- Построение сплайна (вызывается конструктором)
|
||
procedure Spline;
|
||
//
|
||
// Функция 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);
|
||
Spline;
|
||
end;
|
||
|
||
/// Вычисление значения построенного
|
||
/// конструктором сплайна в заданной точке
|
||
function Value(u:real):real;
|
||
//
|
||
// портированная с фортрана функция Seval)
|
||
begin
|
||
var i,j:integer;
|
||
if u<P[0].x then i:=0
|
||
else if u>P[n-1].x then i:=n-2
|
||
else begin
|
||
j:=n-1;
|
||
repeat
|
||
var k:=(i+j) div 2;
|
||
if u<P[k].x then j:=k else i:=k
|
||
until j-i<=1
|
||
end;
|
||
var dx:=u-P[i].x;
|
||
Result:=P[i].y+dx*(aK[0,i]+dx*(aK[1,i]+dx*aK[2,i]))
|
||
end;
|
||
|
||
/// Вычисляет значение первой и второй производной в заданной точке
|
||
function Diff(u:real):(real,real);
|
||
//
|
||
// (модифицированный метод Value)
|
||
begin
|
||
var i,j:integer;
|
||
if u<P[0].x then i:=0
|
||
else if u>P[n-1].x then i:=n-2
|
||
else begin
|
||
j:=n-1;
|
||
repeat
|
||
var k:=(i+j) div 2;
|
||
if u<P[k].x then j:=k else i:=k
|
||
until j-i<=1
|
||
end;
|
||
var dx:=u-P[i].x;
|
||
Result:=(aK[0,i]+dx*(2*aK[1,i]+3*dx*aK[2,i]),2*ak[1,i]+6*dx*aK[2,i])
|
||
end;
|
||
|
||
end;
|
||
{$endregion}
|
||
|
||
{$region Svenn}
|
||
type
|
||
/// Поиск начального интервала неопределенности для локализации
|
||
/// минимума унимодальной функции методом Свенна.
|
||
Svenn=class
|
||
|
||
private
|
||
f:real->real;
|
||
xn,t:real;
|
||
|
||
public
|
||
|
||
constructor(f:real->real; xn,t:real);
|
||
begin
|
||
Self.f:=f;
|
||
Self.xn:=xn;
|
||
Self.t:=t
|
||
end;
|
||
|
||
/// Cобственно, само решение и результат
|
||
function Value:(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) then begin
|
||
if d=t then a:=x2 else b:=x2;
|
||
k+=1; x2:=x1
|
||
end
|
||
until y1>=y2;
|
||
if d=t then Result:=(a,x1,0)
|
||
else Result:=(x1,b,0)
|
||
end;
|
||
|
||
end;
|
||
{$endregion}
|
||
|
||
{$region Zeroin}
|
||
type
|
||
/// Нахождение действительного нуля функции на интервале изоляции
|
||
Zeroin=class
|
||
|
||
private
|
||
tol:real;
|
||
f:real->real;
|
||
|
||
public
|
||
|
||
constructor(f:real->real; tol:real);
|
||
begin
|
||
Self.tol:=tol;
|
||
Self.f:=f;
|
||
end;
|
||
|
||
/// Собственно, само решение и результат
|
||
function Value(ax,bx: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)<Abs(fb) then
|
||
begin (a,b):=(b,c); c:=a; (fa,fb):=(fb,fc); fc:=fa end;
|
||
tol1:=2*eps*Abs(b)+0.5*tol;
|
||
var xm:=0.5*(c-b);
|
||
if (Abs(xm)<tol1) or (fb=0) then begin Result:=b; Exit end;
|
||
if (Abs(e)<tol1) or (Abs(fa)<=Abs(fb)) then (d,e):=(xm,d)
|
||
else begin
|
||
var p,q,s:real;
|
||
if a=c then begin s:=fb/fa; p:=2*xm*s; q:=1-s end
|
||
else begin
|
||
var r:=fb/fc; q:=fa/fc; s:=fb/fa; p:=s*(2*xm*q*(q-r)-(b-a)*(r-1));
|
||
q:=(q-1)*(r-1)*(s-1)
|
||
end;
|
||
if p>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;
|
||
|
||
end;
|
||
{$endregion}
|
||
|
||
//----------------------------------------
|
||
///--
|
||
procedure __InitModule__;
|
||
begin
|
||
end;
|
||
|
||
///--
|
||
procedure __FinalizeModule__;
|
||
begin
|
||
end;
|
||
end. |