Вычислительные методы алгебры и оценивания. И.В. Семушкин (2011) (1185350), страница 48
Текст из файла (страница 48)
е. число неизвестных n;(2) степень переопределенности задачи, т. е. число p, указывающее, восколько раз число уравнений m больше числа неизвестных, m = pn;(3) степень несовместности системы, т. е. вещественное положительноечисло c, показывающее среднеквадратическое значение элементов случайного разностного вектора d = b − Ax̄, где x̄ — нормальное псевдорешение системы Ax ≈ b, относительно среднего (единичного) значения;(4) способ генерации матрицы A.В.
Провести вычислительный эксперимент, используя в нем:(1) десять значений n, n = 1, 2, . . . , 10;(2) три значения p, p = 10, 100 и 1000;(3) три значения c, c = 1/10, 1 и 10;(4) четыре способа генерации матрицы A (см. п. 2 в подразд. 13.14).Результаты эксперимента вывести на экран в виде следующей таблицы:nВычислительный экспериментp=(значение), c=(значение), A=(способ)ПогрешностьПамятьЧисло операцийабсолютная относительная КБайт теоретическое реальноеПри подсчете числа операций учитывать: сложение, умножение, делениеи извлечение квадратного корня. К отчету о работе приложить расчетные30613.14 Задание на лабораторный проект № 9формулы числа операций отдельно по этим видам операций и их сумму.В таблицу выводить только суммарное число операций.Г.
Выполнить отладку программы и продемонстрировать результатыотладки на решении следующей тестовой задачи [82]:hi.Ax ≈ b , A = A(m, 2) = ai1 .. ai2 , i = 1, 2, . . . , m ;m = 4, 8, 12, 16, 20, 24, 28, 32, 36, 40 ;ai1 = wi = sin(2πi/m), ai2 = wi−1 ;bT = (b1, . . . , bm), bi = 2 cos(2πi/m) .Известно (стр. 267), что решением этой задачи является вектор x̄T = (x̄1, x̄2),x̄1 = 2 ctg(2π/m),x̄2 = −2 cosec(2π/m).(13.26)Для демонстрации процесса отладки вывести на экран еще одну таблицу:Отладка программыПогрешностьmабсолютная относительнаяД.
Во всех случаях для оценки абсолютной погрешности использоватьнорму вектора e = x − x̄ видаkek∞ = max |ei | ,iгде x — вычисленное решение, а относительную погрешность определять поформулеδ = kek∞ /kx̄k∞ ,в которой x̄ — точное МHК-решение (нормальное псевдорешение задачи).В отладочной (тестовой) задаче по п. Г решением является вектор (13.26), ав задачах вычислительного эксперимента по п. В — вектор x̄ = (1, 2, . .
. , n).2. Генерация задач для вычислительного экспериментаКак отмечено, для этих задач точное МHК-решение следует задать в видевектора x̄ = (1, 2, . . . , n). Затем следует сгенерировать матрицу A (см. ниже)и образовать вектор b̂ = Ax̄. К нему нужно добавить случайный векторd = cξ, где c — число (см. выше п. 1), а ξ ∼ N (0, 1) — вектор случайныхчисел (от подпрограммы псевдослучайных чисел), взятых из стандартного30713 Устойчивые алгоритмы фильтрациинормального распределения с нулевым средним значением и единичной дисперсией. В результате получаем вектор b = Ax̄ + d для системы уравненийAx ≈ b.Для генерации матрицы A необходимо предусмотреть один из четырехспособов:Способ 1 — матрица A = A(m, n) заполняется случайными числамииз равномерного распределения в диапазоне [−100, +100].
Условно запишемэто так:A = [ Random(m × n) ] .Способ 2 — верхняя часть матрицы A, а именно, ее первые n строк,заполняются по способу 1, а остальные строки образуют подматрицу, в которой только первые q столбцов, где q = [n/2] — ближайшее снизу целое числок n/2, заполняются как в способе 1, а все остальные столбцы — нулевые.Таким образом, матрица имеет вид:A=Random(n × n)Random0Способ 3 — первая часть матрицы A должна формироваться как в способе 2, а остальная часть образуется располагающимися последовательновниз блоками из единичных матриц I размера n × n:Random(n × n)I(n × n)A=···I(n × n)Способ 4 — верхняя часть матрицы A должна формироваться как вспособе 2, остальная же часть строится подобно способу 2, но ненулевые qстолбцов нижней подматрицы заполняются не случайными числами, а располагающимися последовательно вниз блоками из единичных матриц I размера (q × q):Random(n × n)I(q × q)A=···0I(q × q)Замечание 13.18.Hе нужно генерировать всю матрицу A, также как и весь вектор b и весь вектор d, единовременно.
Матрицу A нужно30813.15 Варианты задания на лабораторный проект № 9генерировать построчно, а векторы b и d — поэлементно. Hапример, еслиaT = (a1 , a2, . . . , an ) есть текущая строка матрицы A, то z = aT x̄ + v, гдеz — текущий элемент вектора b, v — текущий элемент вектора d, v = cξ,ξ ∼ N (0, 1) — текущее случайное число из стандартного нормального распределения. Таким образом, последовательно генерируемые данные (aT , z)нужно вводить в алгоритм решения, а также использовать в нем значение r = c2 , имеющее смысл дисперсии ошибки измерения вектора b̂ = Ax̄,исключительно последовательно.13.15Варианты задания на лабораторный проект № 9Общее число вариантов составляет 26 (учитывая подварианты — различия в методах ортогонализации в некоторых из вариантов).Замечание 13.19.
Для всех ковариационных алгоритмов (варианты1–13) в качестве начальных значений можно взять: x0 = 0, P0 = (1/ε2)I,ε → 0.Вариант 1. Стандартный ковариационный алгоритм (Калмана).Найдите его на стр. 274, подразд. 13.3.Вариант 2. Стабилизированный ковариационный алгоритм (Джозефа).Найдите его на стр. 275, подразд. 13.4:(ii) Обработка наблюдений (очередные данные z = aT x̄ + v):α = aT P̃ a + r,K = P̃ a/α,P̂ = (I − KaT )P̃ (I − aK T ) + rKK T .(13.27)Замечание 13.20. Вычислительные затраты существенно зависят отспособа программирования выражений.
Hапример, выражение (13.27) для P̂может быть запрограммировано в следующей последовательности:W1 = I − KaT ,W2 = W1P̃ ,P̂ = W2W1T + r(KK T )(n × n)-матрица(n × n)-матрица30913 Устойчивые алгоритмы фильтрацииили, эквивалентно, в виде:v1 = P̃ a,n-векторP1 = P̃ − Kv1T ,v2 = P1 a,(n × n)-матрицаn-векторP̂ = (P1 − v2 K T ) + (rK)K T ,и в обоих способах можно экономить вычисления, благодаря симметрииматрицы P̂ . Однако второй способ имеет на порядок меньше вычислений:в первом способе выполняется (1, 5n3 + 3n2 + n) умножений, а во второмтолько (4n2 + 2n) умножений.Вариант 3. Квадратно-корневой ковариационный алгоритм Поттера.Найдите его на стр. 276, подразд.
13.5, в следующем виде.(i) Инициализация (начальные значения x0, P0 ):x̃ := x0 ,1/2S̃ := P0 .(ii) Обработка наблюдений (очередные данные z = aT x̄ + v):pf = S̃ T a, α = f T f + r, γ = 1/(1 + r/α) ,K = S̃f /α,Ŝ = S̃ − γKf T ,x̂ = x̃ + K(z − aT x̃) .(iii) Экстраполяция (между повторениями этапа (ii)):S̃ := Ŝ,x̃ := x̂ .Замечание 13.21. Вариант 3, в котором вместо S̃ := Ŝ предусмотрена процедура триангуляризации S̃ := triang Ŝ, дает следующие версии:– Версия 3.1 : обе матрицы, S̃ и Ŝ, — нижние треугольные (S ≡ L), или– Версия 3.2 : обе матрицы, S̃ и Ŝ, — верхние треугольные (S ≡ U ).Именно для этого этап (iii) должен содержать, вместо S̃ := Ŝ, процедурутриангуляризации S̃ := triang Ŝ, матрицы Ŝ. Возможны четыре алгоритмаэтой процедуры: (1) отражения Хаусхолдера, (2) вращения Гивенса, (3) классическая Грама–Шмидта ортогонализация и (4) модифицированная Грама–Шмидта ортогонализация (см.
лабораторный проект № 6). Соответственно31013.15 Варианты задания на лабораторный проект № 9этому, всего имеем 8 подвариантов для указанного варианта 3, сведенных вследующую таблицу:triangS≡L S≡UХаусхолдер 3.1.13.2.1Гивенс3.1.23.2.2ГШО3.1.33.2.3МГШО3.1.43.2.4Вариант 4. Факторизованный LDLT ковариационный алгоритм (Бирмана). Найдите его на стр. 279, подразд. 13.7.Замечание 13.22. Рациональное программирование этого LDLTфакторизованного алгоритма должно экономить память компьютера.
Здесьb поверх De и столбцы Lb поверх столбцов Le (при этомможно записывать Dнужна только поддиагональная часть этих столбцов.Вариант 5. Факторизованный U DU T ковариационный алгоритм (Бирмана). Выведите его самостоятельно аналогично алгоритму варианта 4 [15].Вариант 6. Факторизованный U U T ковариационный алгоритмКарлсона. Выведите его самостоятельно, опираясь на вывод LLT алгоритма, данный в подразд. 13.8. Получите следующий результат [15]:(i) Инициализация (начальные значения x0, P0 ):x̃0 := x0,1/2Ũ := P0 .(ii) Обработка наблюдений (очередные данные z = aT x̂ + v):f = Ũ T a ,α0 = r,f T = (f1, .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
