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

Подпрограмма:  ia20r_c

Назначение

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

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

Для функции 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.

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

    int ia20r_c (integer *n, real *f, integer *ng, real *v1,
            real *v2, integer *mu, real *h, real *p, real *r, real *e,
            integer *isp, integer *lu, real *u, real *fu, integer *kn, real *b)

Параметры

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_c - наилучшая среднеквадратичная аппроксимация одномерной дискретной функции на множестве кусочно - монотонных функций.

Замечания по использованию

  1. 

Кpоме искомого приближения функции  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.

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

int main(void)
{
    /* Initialized data */
    static float v1[11] = { 0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f };
    static float v2[11] = { 1.1f,1.1f,1.1f,1.1f,1.1f,1.1f,1.1f,1.1f,1.1f,
                            1.1f,0.f };
    static int mu[11] = { 1,1,1,1,1,-1,-1,-1,-1,0,0 };
    static float h__[11] = { .1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,.1f,0.f };
    static float p[11] = { 1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f,1.f };
    static float r__[12] = { 0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f };
    static float u[11] = { 0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f };

    /* Builtin functions */
    double sin(double);

    /* Local variables */
    extern int ia20r_c(int *, float *, int *, float *, float *, int *,
                       float *, float *, float *, float *, int *,
                       int *, float *, float *, int *, float *);
    static float b[100], e, f[11];
    static int j, n, i;
    static float s, y, sf, fu;
    static int ng, kn, lu, isp;
    int i__1;
    float r__1;

    n = 11;
    ng = 3;
    lu = 1;
    e = 1e-5f;
    isp = 20;
    u[0] = -.45833333333333331f;
    sf = .025f;
    y = -.5f;
    i__1 = n;
    for (j = 1; j <= i__1; ++j) {
        s = sf * ((float)sin((float)(j * 123.456f + 789.f)) * 2.f);
/* Computing 3rd power */
        r__1 = y;
        f[j - 1] = y - r__1 * (r__1 * r__1) / 3.f + s;
/* l1: */
        y += h__[j - 1];
    }
    ia20r_c(&n, f, &ng, v1, v2, mu, h__, p, r__, &e, &isp, &lu, u, &fu,
            &kn, b);

    printf("\n %16.7e \n",fu);
    for (i = 1; i <= 11; ++i) {
        printf("\n %16.7e \n",u[i-1]);
    }
    return 0;
} /* main */


Результат:

       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)