Документ взят из кэша поисковой машины. Адрес
оригинального документа
: http://num-anal.srcc.msu.su/lib_na/cat/e/ev01r.htm
Дата изменения: Tue Sep 30 16:35:56 2014 Дата индексирования: Sat Apr 9 23:48:51 2016 Кодировка: Windows-1251 |
Текст подпрограммы и версий ( Фортран ) ev01r.zip |
Тексты тестовых примеров ( Фортран ) tev01r.zip |
Текст подпрограммы и версий ( Си ) ev01r_c.zip |
Тексты тестовых примеров ( Си ) tev01r_c.zip |
Текст подпрограммы и версий ( Паскаль ) ev01r_p.zip |
Тексты тестовых примеров ( Паскаль ) tev01r_p.zip |
Решение интегрального уравнения Вольтерра со стационарным ядром методом рягуляризации первого порядка без выбора параметра регуляризации.
Данная подпрограмма реализует алгоритм, изложенный в [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. |
SUBROUTINE EV01R (A, F, N, L, ALFA, H, B, C, U, 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 - | подпрограмма выдачи диагностических сообщений при работе подпрограммы EV01R. |
Замечания по использованию: нет
Решается интегральное уравнение 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.
Параметр регуляризации α = 10 - 8
REAL R(50), F(50), G(100), U(50), WORK(1325), K(50) DATA F /0.025, 0.215, 0.405, 0.595, 0.785, 0.975, 1., 0.92, * 0.71, 0.5, 0.29, 0.08, 0., 0.08, 0.29, 0.5, 0.71, 0.92, * 1., 0.975, 0.785, 0.595, 0.405, 0.215, 0.025, 25*0./ N = 50 ALFA = 1.E - 8 H = 2.E - 2 PI = 3.14159265 N1 = N - 1 DO 11 I = 1, N G(N + I) = 0 S = EXP( - (PI*PI*I*H/8) ) DO 12 J = 1, 15 RAB = EXP( - (2*J + 1)**2*PI*PI*I*H/8 ) 12 S = S + (- 1)**J*(2*J + 1)*RAB 11 R(I) = 0.5*PI*S K(1) = R(1)*H/2 K(N) = R(N)*H/2 DO 13 I = 2, N1 13 K(I) = R(I)*H DO 14 I = 1, N S = 0 DO 15 J = 1, I 15 S = S + F(J)*K(I - J + 1) 14 G(I) = S L = 1 CALL EV01R (R, G, N, L, ALFA, H, WORK(N + 1), WORK, U, IERR) Результаты: 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