Текст подпрограммы и версий
ia20r_p.zip 		Тексты тестовых примеров
tia20r_p.zip

Подпрограмма:  IA20R (модуль IA20R_p)

Назначение

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

Математическое описание

Для функции f (x), a ≤ x ≤ b, строится среднеквадратичное приближение u (x) в классе кусочно - выпуклых функций с ограничениями на первую производную: v1 (x) ≤ u' (x) ≤ v2 (x). Эта задача сводится к минимизации функционала 

                 b  
      min     ∫  p(x)*( u(x) - f(x) )2 dx ,
    v,u(a)    a
                              x
     u(x) = u(a)  +   ∫  v(t) dt
                            a

Hа множестве кусочно - монотонных функций  v (x) = u' (x),  v1 (x) ≤ v (x) ≤ v2 (x), a ≤ x ≤ b. Здесь  p (x) > 0,  v1 (x),  v2 (x) - заданные функции. Значение  u (a) может быть как известным, так и подлежащим определению.

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

                         N
   (1)     min      ∑   pi di ( ui - fi )2
             v,ui     i=1

при условии, что вектоp  v = (v1, v2, ..., vN - 1)  удовлетворяет ограничениям 

   (2)      ui  =  u1 + ( h1*v1 + h2*v2 + ... + hi-1*vi-1 ) ,    i = 2, ..., N
   (3)      v1 i ≤ vi ≤ v2 i  ,                       i = 1, ..., N-1
   (4)      mi ( vi+1 - vi ) ≥ 0  ,                i = 1, ..., N-2

где  di > 0 - коэффициенты квадратурной формулы трапеций на сетке 

 w = { xi : a = x1 <...< xi <...< xN = b } ,   fi = f( xi ) ,
 pi = p( xi ) ,   v1 i = vi( xi ) ,    v2 i = v2( xi ) ,

hi = xi + 1 - xi - шаг сетки  w, параметры  mi определяются условиями кусочной монотонности функции  v = u' и принимают значения  1, 0 или  -1.

Численное решение задачи (1) - (4) основано на методе проекции сопряженных градиентов, который позволяет проводить итерационный процесс минимизации в конечномерном пространстве W1 2, ρ с нормой 

                    N-1
   || v ||   =  [  ∑   pi di vi2  +
                    i=1   
                                         N-1
                                     +  ∑   ρi ( vi - vi-1 )2 ]1/2 ,
                                         i=2

ρi ≥ 0 - заданные числа. В случае достаточной гладкости функции  v (x) задание параметров  ρi > 0 позволяет увеличить точность и улучшить качество приближения для производной. Итерационный процесс считается оконченным после выполнения заданного числа итераций или после достижения заданной точности. 

1. В.А.Моpозов, Н.Л.Гольдман, М.К.Самарин. Метод дескриптивной
   регуляризации и качество приближенных решений,  ИФЖ, T.33,
   N6, 1977.

2. Н.Л.Гольдман. Приближенное решение интегрального уравнения
   Фредгольма I рода в классе кусочно - выпуклых функций
   с ограниченной  первой  производной. Сб. "Численный  анализ на
   ФОРТРАНе, Методы и алгоритмы." Изд-во МГУ, 1978.

3. Ф.П.Васильев. Лекции по методам решения  экстремальных задач.
   Изд-во МГУ, 1974.

Использование 

procedure IA20R(N :Integer; var F :Array of Real; NG :Integer;
                var V1 :Array of Real; var V2 :Array of Real;
                var MU :Array of Integer; var H :Array of Real;
                var P :Array of Real; var R :Array of Real; E :Real;
                ISP :Integer; LU :Integer; var U :Array of Real;
                var FU :Real; var KN :Integer; var B :Array of Real);


   Параметры
N - заданное число значений приближаемой функции  f (тип: целый);
F - заданный вещественный вектоp длины N с компонентами  f ( xi ) ;
NG - заданный признак наличия ограничений на производную (тип: целый):
NG = 0 - если ограничения отсутствуют,
NG = 1 - если выполнены ограничения (3) ,
NG = 2 - если выполнены ограничения (4) ,
NG = 3 - если выполнены ограничения (3) и (4) ;
V1 - вещественный вектоp длины N, первые N - 1 компонент которого содержат нижние ограничения  v1 i ;
V2 - вещественный вектоp длины N, первые N - 1 компонент которого содержат верхние ограничения  v2 i ;
MU - целый вектоp длины N, первые N - 2 компонент которого содержат параметры  mi ;
H - вещественный вектоp длины N, первые N - 1 компонент которого содержат шаги  hi ;
P - вещественный вектоp длины N, содержащий весовые коэффициенты  pi > 0 ;
R - вещественный вектоp длины N + 1 значений параметров  ρi таких, что  R (1) = 0.,  R (J) = ρj ≥ 0,  j = 2, ..., N - 1,  R (N) = 0.,  R (N + 1) = 0.;
E - заданная точность вычислений по градиенту (тип: вещественный);
ISP - заданное максимально допустимое число итераций (тип: целый);
LU - параметр, задаваемый из условия (тип: целый);
LU = 1 - если значение  u (a) известно ,
LU = 0 - если значение  u (a) неизвестно ;
U - вещественный вектоp длины N, задающий произвольное начальное приближение; U (1) = U (a) в случае заданного значения  u (a); в результате работы подпрограммы U содержит искомое приближение   u = (u1, ..., un);
FU - вещественная переменная, содержащая вычисленное значение минимизируемого функционала ;
KN - целая переменная, указывающая причину выхода из подпрограммы:
KN = 0 - если выполнено заданное число итераций ,
KN = 1 - если достигнута точность по градиенту ,
KN = 2 - если достигнута точность по функционалу ,
KN = 3 - если достигнута точность по аргументу ;
B - вещественный рабочий вектоp длины 8N + 2 . 
   Версии:  нет 
   Вызываемые подпрограммы
IA01R - наилучшая среднеквадратичная аппроксимация одномерной дискретной функции на множестве кусочно - монотонных функций. 
   Замечания по использованию
  1. 

Kpоме искомого приближения функции  f (x) подпрограмма позволяет получить приближение для ее производной : компоненты вектоpа  v = (v1, ..., vN - 1) pасположены в B (1), ...,B (N - 1).

  2. 

B случае NG = 0 или NG = 2 массивы V1 и V2 в вычислениях не используются и для них не следует резирвировать память. Соответствующими фактическими параметрами могут быть произвольные идентификаторы. Это замечание относится и к массиву MU в случае NG = 0 или NG = 1.

  3.  Величины, определяющие точность вычисления по аргументу и по функционалу, заданы в самой подпрограмме согласованно с величиной E.

Пример использования 

Unit TIA20R_p;
interface
uses
SysUtils, Math, { Delphi }
Lstruct, Lfunc, UtRes_p, IA20R_p;

function TIA20R: String;

implementation

function TIA20R: String;
var
N,NG,LU,ISP,J,_i,KN :Integer;
E,SF,Y,S,FU :Real;
F :Array [0..10] of Real;
B :Array [0..99] of Real;
const
V1 :Array [0..10] of Real = ( 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 );
V2 :Array [0..10] of Real = ( 1.1,1.1,1.1,1.1,1.1,1.1,1.1,1.1,1.1,1.1,0.0 );
MU :Array [0..10] of Integer = ( 1,1,1,1,1,-1,-1,-1,-1,0,0 );
H :Array [0..10] of Real = ( 0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.1,0.0 );
P :Array [0..10] of Real = ( 1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0,1.0 );
R :Array [0..11] of Real = ( 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 );
U :Array [0..10] of Real = ( 0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0 );
label
_1;
begin
Result := '';  { результат функции }
N := 11;
NG := 3;
LU := 1;
E := 1.E-5;
ISP := 20;
U[0] := -11.0/24.0;
SF := 0.025;
Y := -0.5;
for J:=1 to N do
 begin
  S := SF*(2.0*Sin(123.456*J+789.0));
  F[J-1] := Y-(IntPower(Y,3))/3.0+S;
_1:
  Y := Y+H[J-1];
 end;
IA20R(N,F,NG,V1,V2,MU,H,P,R,E,ISP,LU,U,FU,KN,B);
Result := Result + Format('%s',['   ФYHKЦИOHAЛ=']);
Result := Result + Format('%20.16f ',[FU]) + #$0D#$0A;
Result := Result + Format('%s',['   РЕШЕНИЕ U=' + #$0D#$0A]);
Result := Result + #$0D#$0A;
for _i:=0 to 10 do
 begin
  Result := Result + Format('%20.16f ',[U[_i]]);
  if ( ((_i+1) mod 3)=0 )
   then Result := Result + #$0D#$0A;
 end;
Result := Result + #$0D#$0A;
UtRes('TIA20R',Result);  { вывод результатов в файл TIA20R.res }
exit;
end;

end.

Результат:

       U  =  (-4.58E-1,  -3.97E-1,  -3.00E-1,  -2.04E-1,  -1.07E-1,  4.87E-4,
                 1.08E-1,   2.02E-1,   2.96E-1,   3.91E-1,  4.13E-1)