Текст подпрограммы и версий ia20r_c.zip |
Тексты тестовых примеров tia20r_c.zip |
Наилучшая среднеквадратичная аппроксимация одномерной функции на множестве кусочно - выпуклых функций с ограниченной первой производной.
Для функции 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)