|
Текст подпрограммы и версий ev01r_c.zip |
Тексты тестовых примеров tev01r_c.zip |
Решение интегрального уравнения Вольтерра со стационарным ядром методом рягуляризации первого порядка без выбора параметра регуляризации.
Данная подпрограмма реализyeт aлгоритм, изложенный в [1].
Численное решение интегрального уравнения Вольтерра 1 рода вида
T
(1) ∫ k (t - x) u (x) dx = f (t)
t
t
( ∫ k (t - x) u (x) dx = f (t) )
0
t ∈ [ 0, T ]
на равномерной сетке ti = i H, i = 1, 2, ..., n, H = Т/n, сводится к задаче минимизации
(2) || K u - f || E - min u ,
где K - верхняя (нижняя) треугольная теплицева матрица размера n на n, полученная при дискретизации (1): Ki j = H * k (H ( i - j )).
Используя метод регуляризации А.Н.Тихонова [2], приходим к следующей задаче
(3) || K u - f || E2 + α || D u || E2 - min u ,
где α > 0 - заданный параметр регуляризации, D - дискретный аналог оператора дифференцирования d /dx.
Задача (3) эквивалентна задаче
|| | K | | f | ||
(4) || | | u - | | || - min u ,
|| | α1/2 D | | 0 | || E
для решения которой при помощи метода вращений [3] определим ортогональную матрицу Q, такую, что
| K | | Kα | | f | | f1 |
QT | α1/2 D | = | 0 | , QT | 0 | = | f2 |
где Kα - верхняя (нижняя) треугольная матрица. Тогда решение задачи (4) сводится к решению системы линейных алгебраических уравнений с треугольной матрицей:
(5) Kα u = f1 ,
котоpое осуществляется методом обратной подстановки.
| 1. | L.Elden. An efficient algorithm for the regularization of ill - conditioned least squares problems with triangular toeplitz matrix. - REPORT LITH - MAT - R - 1981 - 25, Sweden, 1981. |
| 2. | А.Н.Тихонов, В.Я.Арсенин. Методы решения некорректных задач. - М., Наука, 1979. |
| 3. | В.В.Воеводин. Численные методы алгебры. - M., Hаука, 1966. |
int ev01r_c(real *a, real *f, integer *n, integer *l,
real *alfa, real *h, real *b, real *c, real *u, integer *ierr)
Параметры
| a - | вещественный вектоp длины n, содержащий первую стpоку (первый столбец) теплицевой матрицы K; |
| f - | вещественный вектоp длины 2n, содержащий в первых n ячейках компоненты заданного вектоpа правой части; последующие n ячеек используются как рабочие; |
| n - | заданная размерность вектоpов f и u и число узлов сетки разбиения (тип: целый); |
| l - | заданный параметр - указатель типа матрицы K: |
|
l = 0, если матрица верхняя треугольная, l = 1, если матрица нижняя треугольная (тип: целый); |
| alfa - | заданный параметр регуляризации (тип: вещественный); |
| h - | заданный шаг разбиения сетки (тип: вещественный); |
| b - | вещественный вектоp длины n (n + 1)/2, используемый как рабочий; |
| c - | вещественный вектоp длины n, используемый как рабочий; |
| u - | вещественный вектоp длины n, содержащий на выходе из подпрограммы вычисленное решение интегрального уравнения в узлах заданной сетки; |
| ierr - | целая переменная, сообщающая об ошибках, обнаруженных в ходе вычислений; |
| ierr= 0 - | если ошибок не было обнаружено; |
| ierr=65 - | в случае, если при заданном значении alfa матрица системы вырождается. |
Версии: нет
Вызываемые подпрограммы
| utev10_c - | подпрограмма выдачи диагностических сообщений при работе подпрограммы ev01r_c. |
Замечания по использованию: нет
Решается интегральное уравнение
t
∫ K (t - s) f (s) ds = g (t) , 0 ≤ t ≤ 1
0
с ядром
∞
K(t) = 0.5 π ∑ (- 1)n (2nh ) exp (-(2n + 1)2 π2 t /8) .
n=0
Точное значение функции f (s) задано таблицей. Используя формулу трапеции, вычисляем приближенное значение правой части g (t), затем по g (t) восстанавливаем значение функции f (s) = u с помощью программы ev01r_c.
Параметр регуляризации α = 10 - 8
int main(void)
{
/* Initialized data */
static float f[50] = { .025f,.215f,.405f,.595f,.785f,.975f,1.f,.92f,.71f,
.5f,.29f,.08f,0.f,.08f,.29f,.5f,.71f,.92f,1.f,.975f,.785f,.595f,
.405f,.215f,.025f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,0.f,
0.f,0.f,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 exp(double);
int pow_ii(int *, int *);
/* Local variables */
static float alfa;
extern int ev01r_c(float *, float *, int *, int *,
float *, float *, float *, float *, float *, int *);
static int ierr;
static float work[1325], g[100], h__;
static int i__, j;
static float k[50];
static int l, n;
static float r__[50], s, u[50];
static int n1;
static float pi, rab;
int i__1, i__2;
n = 50;
alfa = 1e-8f;
h__ = .02f;
pi = 3.14159265f;
n1 = n - 1;
i__1 = n;
for (i__ = 1; i__ <= i__1; ++i__) {
g[n + i__ - 1] = 0.f;
s = (float)exp((float)(-(pi * pi * i__ * h__ / 8)));
for (j = 1; j <= 15; ++j) {
/* Computing 2nd power */
i__2 = (j << 1) + 1;
rab = (float)exp((float)(-(i__2 * i__2) * pi * pi * i__ * h__ / 8));
/* l12: */
s += pow_ii(&c_n1, &j) * ((j << 1) + 1) * rab;
}
/* l11: */
r__[i__ - 1] = pi * .5f * s;
}
k[0] = r__[0] * h__ / 2;
k[n - 1] = r__[n - 1] * h__ / 2;
i__1 = n1;
for (i__ = 2; i__ <= i__1; ++i__) {
/* l13: */
k[i__ - 1] = r__[i__ - 1] * h__;
}
i__1 = n;
for (i__ = 1; i__ <= i__1; ++i__) {
s = 0.f;
i__2 = i__;
for (j = 1; j <= i__2; ++j) {
/* l15: */
s += f[j - 1] * k[i__ - j];
}
/* l14: */
g[i__ - 1] = s;
}
l = 1;
ev01r_c(r__, g, &n, &l, &alfa, &h__, &work[n], work, u, &ierr);
i__1 = n;
for (i__ = 1; i__ <= 21; i__ += 20) {
printf("\n %16.7e %16.7e %16.7e %16.7e ",
u[i__-1], u[i__+4], u[i__+9], u[i__+14]);
}
printf("\n %16.7e %16.7e \n",u[40], u[45]);
printf("\n %16.7e \n",alfa);
return 0;
} /* main */
Результаты:
alfa = 1.00000000-08
| 4.18939426-02 9.47262731-01 2.69209987-01 4.94652786-01
u = | 8.07735236-01 9.34217208-04 8.94488047-03 -2.22705150-03
| -3.81489539-04 5.87342988-03