Created
April 13, 2013 06:32
-
-
Save ckazzku/5377323 to your computer and use it in GitHub Desktop.
Решение уравнения Фредгольма I рода методом итеративной регуляризации Фридмана с выбором числа итераций способом обобщенной невязки с использованием значений погрешностей правой части и оператора
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # -*- coding: utf-8 -*- | |
| from __future__ import division | |
| from math import sqrt, exp, pi | |
| def friedman(f, node_x, node_s, delta_f, ksi, sigma, max_iterations, kernel): | |
| """ | |
| Решение уравнения Фредгольма I рода методом итеративной регуляризации | |
| Фридмана с выбором числа итераций способом обобщенной невязки | |
| с использованием значений погрешностей правой части и оператора | |
| Входные параметры: | |
| f - правая часть (массив длиной M) | |
| node_x - узлы по координате х (массив длиной M) | |
| node_s - узлы по координате s (массив длиной N) | |
| delta_f - поточечная среднеквадратическая погрешность правой части | |
| ksi - погрешность оператора | |
| sigma - параметр метода, 0 < sigma < 2 | |
| max_iterations - максимальное число итераций | |
| kernel - ядро интеграла | |
| Рабочие массивы: | |
| p - массив длиной M | |
| r - массив длиной N | |
| z - двухмерный массив длиной M x N | |
| g - массив длиной N(N + 1)/2 | |
| ff - массив длиной N | |
| ym - массив длиной N | |
| ym1 - массив длиной N | |
| conditions_unsatisfied - индикатор ошибки: | |
| 0 - выполнены все условия | |
| 1 - нарушено хотя бы одно из условий (E1)-(E4) | |
| 2 - нарушено хотя бы одно из условий (E5)-(E6) | |
| 3 - объединение случаев 1 и 2 | |
| (E1): N >= N_min | |
| (E2): M >= M_min | |
| (E3): монотонность сетки узлов по x | |
| (E4): монотонность сетки узлов по s | |
| (E5): delta_f >= 0, ksi >= 0 | |
| (E6): 0 < sigma < 2, max_iterations > 1 | |
| """ | |
| len_x, len_s = len(node_x), len(node_s) | |
| g = [0 for i in range(len_s * (len_s + 1) // 2)] | |
| ym1 = [0 for i in range(len_s)] | |
| ym = [0 for i in range(len_s)] | |
| y = [0 for i in range(len_s)] | |
| conditions_unsatisfied = not (is_monotonic(node_x) and \ | |
| is_monotonic(node_s)) and 1 or 0 | |
| if not (delta_f > 0 and ksi >= 0 and 0 < sigma < 2 and \ | |
| max_iterations >= 1): | |
| conditions_unsatisfied += 2 | |
| if conditions_unsatisfied: | |
| result = [None for i in range(7)] | |
| else: | |
| result = [0 for i in range(7)] | |
| delta = delta_f * sqrt(node_x[len_x - 1] - node_x[0]) | |
| # Вычисление p(i), r(j), z(i,j) = r(j)K(i,j), G(k,j), F(k) | |
| p, r = coef(node_x), coef(node_s) | |
| z = eval_z(node_x, node_s, r, kernel) | |
| rho = eval_rho(r) | |
| g = eval_g(p, z, rho) | |
| ff = eval_f(f, p, z, rho) | |
| nu = sigma / eval_lambda(g, len_s) | |
| # оператор с погрешностью или нет | |
| operator_have_imprecision = (ksi != 0) | |
| # 1 этап: y_m, beta_m, norm_ym, mu | |
| mu = discrepancy(f, p, ym, z) | |
| norm_ym_prev = 0 | |
| norm_ym_next = operator_have_imprecision and norm_in_L2(ym, r) or 0 | |
| # Oпределениe mu | |
| for m in range(max_iterations): | |
| ym = solution(g, ff, ym1, nu) | |
| # if m in [0,9,99,999]: | |
| # print 'm=%s\nym=%s' % (m+1,ym) | |
| beta_m = discrepancy(f, p, ym, z) | |
| if operator_have_imprecision: | |
| norm_ym_previous = norm_in_L2(ym, r) | |
| if not (mu < beta_m or operator_have_imprecision and \ | |
| norm_ym_prev < norm_ym_next): | |
| mu = beta_m | |
| if operator_have_imprecision: | |
| norm_ym_next = norm_ym_prev | |
| ym1 = [ym[i] for i in range(len_s)] | |
| delta2_plus_mu = delta**2 + mu | |
| norm_f = norm_in_L2(f, p) | |
| # условие существования и единственности оптимального пар-ра рег-ии | |
| optimal_parameter_exist_and_unique = (norm_f**2 > delta2_plus_mu) | |
| # 2 этап: m_d, y_m_d | |
| m = optimal_parameter = 0 | |
| if optimal_parameter_exist_and_unique: | |
| for i in range(len_s): | |
| ym[i] = ym1[i] = 0 | |
| while m < max_iterations: | |
| ym = solution(g, ff, ym1, nu) | |
| beta = discrepancy(f, p, ym, z) | |
| if operator_have_imprecision: | |
| dzeta = (delta + ksi * norm_in_L2(ym, r))**2 + mu | |
| else: | |
| dzeta = delta2_plus_mu | |
| kappa = beta - dzeta | |
| if kappa > 0: | |
| ym1 = [ym[i] for i in range(len_s)] | |
| else: | |
| y = [ym[i] for i in range(len_s)] | |
| result[6] = norm_difference_in_L2(y, ym1, r) | |
| optimal_parameter = m + 1 | |
| m = max_iterations | |
| m += 1 | |
| else: | |
| print 'Optimal parameter not exist or not unique!' | |
| result[6] = 0 | |
| # сбор результатов | |
| norm_y = norm_in_L2(y, r) | |
| result[0] = optimal_parameter | |
| result[1] = sqrt(beta) | |
| if optimal_parameter_exist_and_unique: | |
| beta = discrepancy(f, p, y, z) | |
| result[2] = result[1] / norm_f | |
| else: | |
| beta = norm_f**2 | |
| result[2] = 1 | |
| result[3] = beta - ((delta + ksi * norm_y)**2 + mu) | |
| result[4] = mu | |
| result[5] = norm_y | |
| return (conditions_unsatisfied, result, y) | |
| def is_monotonic(node_x): | |
| """ | |
| Проверка монотонноcти сетки узлов | |
| при задании узлов x(i) неравномерной сетки | |
| """ | |
| len_x = len(node_x) | |
| if len_x < 2: | |
| return False | |
| else: | |
| for i in range(len_x - 1): | |
| if node_x[i] >= node_x[i + 1]: | |
| return False | |
| return True | |
| def coef(node_s): | |
| """ | |
| Вычисление коэффициентов квадратурной формулы трапеций по значениям | |
| узлов в случае неравномерной сетки узлов | |
| """ | |
| len_s = len(node_s) | |
| result = [0 for i in range(len_s)] | |
| result[0] = 0.5 * (node_s[1] - node_s[0]) | |
| if len_s > 2: | |
| for i in range(1, len_s - 1): | |
| result[i] = 0.5 * (node_s[i + 1] - node_s[i - 1]) | |
| result[len_s - 1] = 0.5 * (node_s[len_s - 1] - node_s[len_s - 2]) | |
| return result | |
| def eval_g(p, z, rho): | |
| """ | |
| Вычисление коэффициентов G(k,j) | |
| т.к. матрица семмитричная, то используется только | |
| ее верхний треугольник в виде одномерного массива | |
| """ | |
| len_x, len_s = len(p), len(z[0]) | |
| result = [0 for i in range(len_s * (len_s + 1) // 2)] | |
| for k in range(len_s): | |
| for j in range(k + 1): | |
| result[k * (k + 1) // 2 + j] = sum(p[i] / rho * z[i][k] * z[i][j] \ | |
| for i in range(len_x)) | |
| return result | |
| def eval_f(f, p, z, rho): | |
| """ | |
| Вычисление коэффициентов F(k) | |
| """ | |
| len_x, len_s = len(f), len(z[0]) | |
| result = [sum(p[i] / rho * z[i][j] * f[i] for i in range(len_x)) \ | |
| for j in range(len_s)] | |
| return result | |
| def eval_z(node_x, node_s, r, kernel): | |
| """ | |
| Вычисление коэффициентов z(i,j) = r(j)K(i,j) | |
| """ | |
| len_x, len_s = len(node_x), len(node_s) | |
| result = [[0 for j in range(len_s)] for i in range(len_x)] | |
| for j in range(len_s): | |
| for i in range(len_x): | |
| result[i][j] = r[j] * kernel(node_x[i], node_s[j]) | |
| return result | |
| def eval_rho(r): | |
| """ | |
| Вычисление нормирующего делителя | |
| """ | |
| len_r = len(r) | |
| result = 2 * (r[0] + r[len_r - 1]) | |
| if len_r > 2: | |
| result = sum(r[i] for i in range(1, len_r - 1)) | |
| result /= len_r | |
| return result | |
| def eval_lambda(g, n): | |
| """ | |
| Вычисление максимального собственного значения lambda | |
| положительно определенной симметричной матрицы G | |
| """ | |
| result = tmp = 0 | |
| for i in range(n): | |
| for j in range(n): | |
| tmp = g[i * (i + 1) // 2 + j if i >= j else i + j * (j + 1) // 2] | |
| result += tmp**2 | |
| result = sqrt(result) | |
| return result | |
| def solution(g, ff, known_y, nu): | |
| """ | |
| Вычисление решения y(m) по известному y(m-1) | |
| """ | |
| result = [0 for i in range(len(known_y))] | |
| for i in range(len(ff)): | |
| tmp = ff[i] | |
| for j in range(len(known_y)): | |
| tmp -= g[i * (i + 1) // 2 + j if i >= j \ | |
| else i + j * (j + 1) // 2] * known_y[j] | |
| result[i] = known_y[i] + nu * tmp | |
| return result | |
| def discrepancy(f, p, y, z): | |
| """ | |
| Вычисление невязки beta = ||Ay — f||^2 , используемой в способах | |
| (обобщенной) невязки выбора параметра регуляризации alfa | |
| в методе выбора числа итераций ( beta = beta(m), y = y(m) ) | |
| """ | |
| len_y = len(y) | |
| result = 0 | |
| for i in range(len(f)): | |
| tmp = sum(z[i][j] * y[j] for j in range(len_y)) | |
| result += p[i] * (tmp - f[i])**2 | |
| return result | |
| def norm_in_L2(y, r): | |
| """ | |
| Вычисление нормы вектора y(m) в пространстве L2 | |
| """ | |
| result = sqrt(sum(r[i] * y[i]**2 for i in range(len(r)))) | |
| return result | |
| def norm_difference_in_L2(y1, y2, r): | |
| """ | |
| Вычисление нормы разности векторов y1, y2 в пространстве L2 | |
| """ | |
| result = sqrt(sum(r[i] * (y1[i] - y2[i])**2 for i in range(len(r)))) | |
| return result | |
| def test(): | |
| """ | |
| Тестирование функции friedman ( данные взяты из учебника ) | |
| """ | |
| f = [0.0138, 0.0317, 0.0539, 0.0863, 0.1338, 0.1873, 0.2641, 0.3490,\ | |
| 0.4367, 0.5220, 0.6180, 0.6874, 0.7417, 0.7805, 0.7993, 0.7835,\ | |
| 0.7445, 0.6814, 0.6208, 0.5331, 0.4312, 0.3446, 0.2708, 0.1868,\ | |
| 0.1327, 0.0986, 0.0632, 0.0317, 0.0242] | |
| x = [- 1.4 + 0.1 * i for i in range(29)] | |
| s = [- 1 + 0.2 * i for i in range(11)] | |
| kernel = lambda x, s: sqrt(4 / pi) * exp(-4 * (x - s)**2) | |
| result = friedman(f, x, s, 0.5 * 10**-2, 10**-2, 0.5, 1000, kernel) | |
| print 'Results:\n ier: %s\n ymd params: %s\n ymd: %s' % result | |
| if __name__ == '__main__': | |
| # kernel = lambda x, s: exp(-(x - s)**2 * v**2 / 2) / (sqrt(2 * pi) * v) | |
| test() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment