Вариационные методы решения систем линейных уравнений
Решение дифференциальных уравнений в частных производных. Метод минимальных невязок, минимальных поправок, скорейшего спуска, сопряженных градиентов. Алгоритмы и блок-схемы решения. Руководство пользователя программы. Решение системы с матрицей.
| Рубрика | Математика |
| Вид | курсовая работа |
| Язык | русский |
| Дата добавления | 21.01.2014 |
| Размер файла | 380,3 K |
Отправить свою хорошую работу в базу знаний просто. Используйте форму, расположенную ниже
Студенты, аспиранты, молодые ученые, использующие базу знаний в своей учебе и работе, будут вам очень благодарны.
Размещено на http://www.allbest.ru/
Вариационные методы решения систем линейных уравнений
1. Постановка задачи
Большое количество задач математики и физики сводится к решению дифференциальных уравнений в частных производных, решение которых, в свою очередь, приводит к решению систем линейных алгебраических уравнений (СЛАУ). Методы решения СЛАУ, можно разделить на прямые методы, позволяющие получить точное решение системы (к ним относится, например, метод Гаусса, метод Крамера и т.д.) и итерационные методы, позволяющие получить решение системы с заданным приближением. Среди итерационных методов особое место занимают методы вариационного типа, особенностью которых является то, что при их использовании нет необходимости знать границы спектра матрицы системы.
Итак, сформулируем задачу:
Дана система линейных уравнений
где A симметричная положительно определенная матрица.
Требуется решить данную СЛАУ следующими итерационными методами вариационного вида:
1. Метод минимальных невязок
2. Метод минимальных поправок
3. Метод скорейшего спуска
4. Метод сопряженных градиентов
Методы решения СЛАУ рассматриваются практически в каждой книге по численным методам. Общее описание алгоритмов для каждого метода дано, например, в книге Самарского А.А., Гулина А.В. «Численные методы» [3]. Также популярными книгами являются: М.Ю. Баландин, Э.П. Шурина «Методы решения СЛАУ большой размерности» [1], Y. Saad «Iterative Method for Sparse Linear System» [5]. В перечисленных книгах, описаны такие методы как CG (Conjugate Gradient - метод сопряжённых градиентов) и MinRES (Minimal Residual - метод минимальных невязок). В книге [5] изложены не только базовые алгоритмы данных методов, но и указания к их практической реализации.
2. Описание методов решения задачи
2.1 Метод минимальных невязок
Рассмотрим систему с матрицей A=AT>0. Обозначим через
(1)
невязку, которая получается при подстановке приближенного значения xk, полученного на k-й итерации, в уравнение (1). Заметим, что погрешность zk=xkх и невязка rk связаны равенством Azk=rk.
Рассмотрим явный итерационный метод
, (2) <=> (3)
где параметр k+l выбирается из условия минимума погрешности ||rk+1|| при заданной норме ||rk||. Метод (2) называется методом минимальных невязок.
Найдем явное выражение для параметра k+l. Из (3) получаем
, (4) => (5)
т.е. rk удовлетворяет тому же уравнению, что и погрешность zk=xkх. Возводя обе части уравнения (5) для невязки скалярно в квадрат, получаем
(6)
Отсюда видно, что |||rk+1|| достигает минимума, если
(7)
Таким образом, в методе минимальных невязок переход от k-й итерации к (k+1) - й осуществляется следующим образом. По найденному значению xk вычисляется вектор невязки rk=Axkf и по формуле (7) находится параметр k+l. Затем по формуле (3) досчитывается вектор xk+1.
Метод минимальных невязок (3), (7) сходится с той же скоростью, что и метод простой итерации с оптимальным параметром .
2.2 Метод минимальных поправок
Рассмотрим неявный итерационный метод вида
, (8)
Запишем этот метод в виде
. (9)
Вектор k =B1rk называется поправкой на (k+1) - й итерации. Она удовлетворяет тому же уравнению, что и погрешность zk=xkx неявного метода, т.е. уравнению
. (10)
Будем предполагать, что B=BT>0. Методом минимальных поправок называется неявный итерационный метод (8), в котором параметр k+l выбирается из условия минимума нормы ||k+1||B=(Bk+1,k+1)1/2 при заданном векторе k. В случае B=Е метод минимальных поправок совпадает с методом минимальных невязок.
Найдем выражение для итерационного параметра k+l. Перепишем (10) в виде
и вычислим
.
Отсюда следует, что минимум этого выражения достигается при
. (11)
Для реализации метода минимальных поправок требуется на каждой итерации решить две системы уравнений Bk= rk и Bvk=Ak, откуда найдем вектор vk = B1 Аk (воспользуемся формулой (9)).
Скорость сходимости метода минимальных поправок определяется границами спектра обобщенной задачи на собственные значения
. (12)
2.3 Метод скорейшего спуска
Рассмотрим явный метод (2) и выберем итерационный параметр k+l из условия минимума ||zk+1||A при заданном векторе zk, где zk+1=xk+1x. Поскольку погрешность zk удовлетворяет уравнению
, то
.
Следовательно, это выражение достигает минимума, если
.
Величина zk=xkx здесь неизвестна (так как неизвестно точное решение х), поэтому надо учесть, что Azk=rk=Axkf, и вычислять k+l по формуле
. (13)
2.4. Метод сопряженных градиентов
Пусть A - матрица системы . Рассмотрим следующий класс неявных двухшаговых итерационных методов:
. (14)
Здесь , k+l, k+l - итерационные параметры, которые будут определены далее. Начальное приближение x0 будем задавать произвольно, а вектор x1 вычислять по одношаговой формуле, которая получается из (14) при k=0 и l=1,
. (15)
Если параметры k+l, k+l найдены, то новое приближение xk+l выражается через два предыдущих значения хk и xk1 по формуле
. (16)
Для погрешности zk=xkx из (16) получаем уравнения
.
Введем, как и ранее, вспомогательную функцию vk=A1/2zk, для которой ||vk||=||zk||A. Функция vk удовлетворяет уравнениям
.(17)
где C=A1/2B1A1/2. Будем считать, что A и B - симметричные положительно определенные матрицы, удовлетворяющие неравенствам
. (18)
Исключая последовательно векторы vl, v2, , vkl из уравнений (17) находим, что
, (19)
где Pk(C) - многочлен степени k от оператора С, удовлетворяющий условию Pk(0)=E.
Поставим задачу выбрать итерационные параметры k, k так, чтобы при любом n=1, 2,… была бы минимальной ||vn||=||zn||A.
Параметр 1 находится из условия минимума ||v1|| при заданном векторе v0. Так же, как и в методе скорейшего спуска, получаем
. (20)
Отметим, что при таком выборе 1 выполняется равенство (Cv1, v0)=0, т.е. векторы v1 и v0 ортогональны в смысле скалярного произведения
.
Далее, рассмотрим погрешность (19) и запишем многочлен Pk(C) в виде
, (21)
где ai(k) - числовые коэффициенты, определяемые параметрами i, i, i=l, 2,…, k. Тогда
. (22)
Найдем условия, которым должны удовлетворять коэффициенты ai(k), минимизирующие ||vn||2. Из (22) получаем
, (23)
т.е. ||vn||2 является многочленом второй степени по переменным a1(n),…, an(n). Приравнивая к нулю частные производные ||vn||2/aj(n), j=1, 2,…, n, приходим к системе уравнений
, (24)
решение которой ai(n), i=1, 2,…, n и обращает в минимум ||vn||2.
3. Алгоритмы и блок-схемы решения
3.1 Метод минимальных невязок
1. Задаем вектор х0 (начальное приближение).
2. Положим xk = x0, k = 0 (номер итерации)
3. Вычисляем вектор rk = Axk - b (невязка начального приближения).
4. Вычисляем скаляр ?k+1 = (rk, Ark) / ||Ark||2.
5. Вычисляем вектор xk+1 = xk + ?k+1rk (очередное приближение).
6. Вычисляем вектор rk+1 = rk - ?k+1Ark. (невязка k+1 приближения).
7. Проверяем выполнение неравенства ||r(k+1)||2 =< 0.001; если «да», то останавливаем работу алгоритма и выводим результаты.
8. Положим k = k + 1 и вернемся к шагу 4.
Рис. 1. Блок-схема метода минимальных невязок
3.2 Метод минимальных поправок
1. Задаем вектор х0 (начальное приближение), матрицу B = BT > 0.
2. Положим xk = x0, k = 0 (номер итерации)
3. Вычисляем вектор rk = Axk - b (невязка начального приближения).
4. Вычисляем вектор ?k = rkBk-1.
5. Вычисляем скаляр ?k+1 = (A?k, ?k) / (B-1A?k, A?k) (шаговый множитель).
6. Вычисляем вектор xk+1 = xk - ?k+1 B-1A?k (очередное приближение).
7. Вычисляем вектор rk+1 = Axk+1 - b. (невязка k+1 приближения).
8. Положим k = k + 1.
9. Проверяем выполнение неравенства ||rk||2 =< 0.001; если «да», то останавливаем работу алгоритма и выводим результаты, иначе к шагу 3.
Рис. 2. Блок-схема метода минимальных поправок
3.3 Метод скорейшего спуска
1. Задаем вектор х(0) (начальное приближение).
2. Положим xk = x0, k = 0 (номер итерации)
3. Вычисляем вектор rk = Axk - b (невязка начального приближения).
4. Вычисляем скаляр ?k+1 = (rk, rk) / (Ark, Ark) (шаговый множитель).
5. Вычисляем вектор xk+1 = xk - ?k+1rk (очередное приближение).
6. Вычисляем вектор rk+1 = Axk+1 - b. (невязка k+1 приближения).
7. Положим k = k + 1.
8. Проверяем выполнение неравенства ||rk||2 =< 0.001; если «да», то останавливаем работу алгоритма и выводим результаты, иначе возвращаемся к шагу 3.
Рис. 3. Блок-схема метода скорейшего спуска
3.4 Метод сопряженных градиентов
1. Задаем вектор х(0) (начальное приближение).
2. Положим xk = x0, k = 0 (номер итерации)
3. Вычисляем вектор rk = Axk - b (невязка начального приближения).
4. Положим вектор pk= rk (сопряженный вектор)
5. Вычисляем скаляр ?k = (rk, rk) / (Apk, pk).
6. Вычисляем вектор xk+1 = xk + ?k pk (очередное приближение)
7. Вычисляем вектор rk+1 = rk - ?k Apk. (невязка k+1 приближения)
8. Вычисляем скаляр ?k = (rk+1, rk+1) / (rk, rk).
9. Проверяем выполнение неравенства ||rk+1||2 =< 0.001; если «да», то останавливаем работу алгоритма и выводим результаты
10. Вычисляем pk+1= rk+1+ ?kpk
11. Положим k = k + 1.
12. Возвращаемся к шагу 5.
4. Руководство пользователя
В окне программы пользователю предлагается ввести размерность матрицы системы и требуемую точность, также при желании можно задать величину диагонального преобладания, либо убрать преобладание вообще (см. рис. 5).
Рис. 5. Окно программы
При нажатии на кнопку ОК генерируется симметричная матрица заданной размерности, а также правая часть. Далее, при нажатии на кнопку с названием метода, в соответствующий столбец выводится решение и количество итераций (см. рис. 6).
Рис. 6. Решение системы с матрицей размерностью 10х10
Матрицу можно модифицировать, прибавив к главной диагонали заданную величину, для этого нужно ввести необходимую величину в поле «Преобладание», поставить флаг «Модифицировать» и нажать ОК.
Также существует возможность записи сгенерированной матрицы в файл и чтения существующей матрицы из файла. Пользователь может выбрать для записи и чтения один из трёх файлов.
5. Тестирование и оптимизация
Изначально программа тестировалась на симметричных матрицах с диагональным преобладанием, то есть при построении матрицы на главную диагональ прибавлялась 1000, что давало очень хорошую сходимость при заданной точности 0,001 у всех четырёх методов. При этом было видно, что наихудшую скорость сходимости дает метод минимальных поправок, наилучшую - метод сопряжённых градиентов. Разница в скорости сходимости становилась заметна, при размерности матрицы больше чем 200х200 элементов, что видно на графике (см. рис. 7).
Рис. 7. Зависимость количества итераций от размерности матрицы
Для более наглядной демонстрации работы вариационных методов в программу была добавлена возможность генерировать матрицы без диагонального преобладания, а также с диагональным преобладанием, заданным пользователем. В результате были сделаны следующие наблюдения: при отсутствии диагонального преобладания метод минимальных невязок очень медленно сходится (за несколько сотен тысяч итераций), при этом он сходится только для точности не меньше 0,01 и матриц размерностью не больше 50. Аналогичная ситуация наблюдается с методами минимальных поправок и скорейшего спуска. Если задавать небольшое диагональное преобладание, то скорость сходимости увеличивается. На рисунке 8 показано, как изменяется количество итераций в зависимости от диагонального преобладания при размерности матрицы 10х10 и точности, равной 0,01.
Рис. 8 Зависимость количества итераций от диагонального преобладания
Что касается метода сопряжённых градиентов, то диагональное преобладание незначительно влияет на скорость его сходимости, так, в случае матрицы 10х10 элементов без диагонального преобладания и при точности, равной 0,01, метод минимальных невязок сошёлся за 25329 итераций, а метод сопряжённых градиентов - за 11 итераций.
Рис. 9. Зависимость количества итераций в методе сопряженных градиентов от диагонального преобладания
Таким образом, из всех исследованных методов метод сопряжённых градиентов является наиболее эффективным.
Для хранения матриц и векторов были использованы динамические массивы, что ускоряет работу программы и позволяет пользователю практически неограниченно задавать размерность матрицы системы. Также в методе минимальных поправок матрично-матричное умножение заменено матрично-векторным, что позволяет значительно снизить время, затрачиваемое на вычисления.
Выводы
В ходе выполнения данной курсовой работы были исследованы такие методы решения СЛАУ, как метод минимальных невязок, метод минимальных поправок, метод скорейшего спуска и метод сопряженных градиентов. Составлены алгоритмы решения СЛАУ перечисленными методами, которые в дальнейшем были использованы для написания программы. Сравнение решений, полученных с помощью данных методов, с точным решением показало состоятельность методов. В ходе исследования сходимости в случае матриц с диагональным преобладанием, было обнаружено, что скорость сходимости метода сопряжённых градиентов наибольшая. В случае матриц без диагонального преобладания метод сопряжённых градиентов сходится значительно быстрее всех остальных методов, причём разница в сходимости может составлять десятки тысяч раз. Методы минимальных невязок, минимальных поправок и скорейшего спуска применимы только в случае матриц с хорошим диагональным преобладанием. В отличие от этих методов, метод сопряжённых градиентов сходится для любых симметричных матриц, при этом он даёт настолько лучшую сходимость по сравнению с остальными методами, что это делает его самым эффективным среди рассмотренных методов.
Список источников
1. Баландин М.Ю., Шурина Э.П. «Методы решения СЛАУ большой размерности» - Новосибирск: Изд-во НГТУ, 2000
2. Бахвалов Н.С., Жидков Н.М. «Численные методы» - М.: Наука, 1986
3. Самарский А.А., Гулин А.В. «Численные методы», Изд. Наука, 1989
4. Фаддеев В.К., Фаддеева В.Н. «Вычислительные методы линейной алгебры» - М.: Наука 1958.
5. Saad Y. Iterative methods for sparse linear systems. Boston: PWS Publ. Co., 2000.
Приложение
// -
#include <vcl.h>
#include <time.h>
#include <math.h>
#include <stdio.h>
#pragma hdrstop
#include «Unit1.h»
// -
#pragma package (smart_init)
#pragma resource «*.dfm»
TForm1 *Form1;
int n, m;
float epsilon; // Точность
double**a; // Матрица
double *f; // Столбец правых частей
double *x; // Решение
double *ax;
double**b;
bool ch=0;
double *tm;
// -
__fastcall TForm1:TForm1 (TComponent* Owner)
: TForm(Owner)
{
sg[0]=StringGrid2;
sg[1]=StringGrid3;
sg[2]=StringGrid4;
sg[3]=StringGrid5;
sg[4]=StringGrid6;
sg[5]=StringGrid7;
sg[6]=StringGrid8;
}
// -
// Матрично-векторное умножение
void mv_mult (double **a, double *x, double *ax)
{
for (int i=0; i<n; i++)
{
ax[i]=0;
for (int j=0; j<n; j++)
{
ax[i]+=a[i] [j]*x[j];
}
}
}
// -
// Векторное вычитание
void vv_substr (double *ax, double *f, double *r)
{
for (int i=0; i<n; i++)
{
r[i]=ax[i] - f[i];
}
}
// -
// Векторное сложение
void vv_addit (double *ax, double *f, double *r)
{
for (int i=0; i<n; i++)
{
r[i]=ax[i]+f[i];
}
}
// -
// Норма вектора
double norma (double * vec)
{
double nr=0;
for (int i=0; i<n; i++)
{
nr+=vec[i]*vec[i];
}
return (sqrt(nr));
}
// -
// Скалярное умножение векторов
double vv_mult (double *a, double *x)
{
double ax=0;
for (int i=0; i<n; i++)
{
ax+=a[i]*x[i];
}
return(ax);
}
// -
// Проверка на заданную точность нормы
bool check (double *a)
{
double mod;
mod=norma(a);
if (mod<epsilon) return(1);
else return(0);
}
// -
// Построение матрицы B для метода минимальных поправок
void mkB()
{
ch=1;
b=new double* [n];
for (int i=0; i<n; i++)
{
b[i]=new double[n];
for (int j=0; j<n; j++)
if (i!=j) b[i] [j]=0;
else b[i] [j]=a[i] [j] - 300;
}
}
// -
// Заполнение матрицы
void __fastcall TForm1: Button5Click (TObject *Sender)
{
Button1->Enabled=true;
Button2->Enabled=true;
Button3->Enabled=true;
Button4->Enabled=true;
CheckBox2->Enabled=true;
time_t t;
time(&t); srand((unsigned) t);
n=StrToInt (Edit1->Text);
m=StrToInt (Edit3->Text);
epsilon=StrToFloat (Edit2->Text);
StringGrid1->ColCount=n;
StringGrid1->RowCount=n;
for (int i=0; i<6; i++)
sg[i]->RowCount=n;
f=new double [n];
if (CheckBox2->Checked==false) {
a=new double* [n];
tm=new double [n];
x=new double [n];
for (int i=0; i<n; i++)
a[i]=new double [n];
for (int i=0; i<n; i++)
{
double b=rand()%200-100;
double c=rand()%100+1;
x[i]=b/c;
sg[1]->Cells[0] [i]=FloatToStr (x[i]);
for (int j=0; j<=i; j++)
{
double b=rand()%200-100;
double c=rand()%100+1;
a[i] [j]=b/c;
a[j] [i]=a[i] [j];
if ((i==j)&&(CheckBox1->Checked==true)) a[i] [j]+=m;
if (i==j) tm[i]=a[i] [j];
StringGrid1->Cells[j] [i]=FloatToStr (a[i] [j]);
StringGrid1->Cells[i] [j]=FloatToStr (a[j] [i]);
}
}}
else
for (int i=0; i<n; i++)
{
a[i] [i]=tm[i]+(double) m;
tm[i]+=m;
StringGrid1->Cells[i] [i]=FloatToStr (a[i] [i]);
StringGrid1->Cells[i] [i]=FloatToStr (a[i] [i]);
}
mv_mult (a, x, f);
for (int i=0; i<n; i++)
sg[0]->Cells[0] [i]=FloatToStr (f[i]);
}
матрица вариационный уравнение линейный
Размещено на Allbest.ru
Подобные документы
Решение эллиптических и параболических дифференциальных уравнений в частных производных. Суть метода Кранка-Николсона и теории разностных схем для теплопроводности. Построение численных методов с помощью вариационных принципов, описание Matlab и Mathcad.
курсовая работа [1,4 M], добавлен 13.03.2011Методы решения систем линейных алгебраических уравнений, их характеристика и отличительные черты, особенности и сферы применения. Структура метода ортогонализации и метода сопряженных градиентов, их разновидности и условия, этапы практической реализации.
курсовая работа [197,8 K], добавлен 01.10.2009Метод Гаусса–Жордана: определение типа системы, запись общего решения и базиса. Выражение свободных переменных с использованием матричного исчисления. Нахождение координат вектора в базисе. Решение системы уравнений по правилу Крамера и обратной матрицей.
контрольная работа [200,4 K], добавлен 17.12.2010Основные правила решения системы заданных уравнений методом Гаусса с минимизацией невязки и методом простых итераций. Понятие исходной матрицы; нахождение определителя для матрицы коэффициентов. Пример составления блок-схемы метода минимизации невязок.
лабораторная работа [264,1 K], добавлен 24.09.2014Анализ методов решения систем дифференциальных уравнений, которыми можно описать поведение материальных точек в силовом поле, законы химической кинетики, уравнения электрических цепей. Этапы решения задачи Коши для системы дифференциальных уравнений.
курсовая работа [791,0 K], добавлен 12.06.2010Проверка совместности системы уравнений, ее решение матричным методом. Координаты вектора в четырехмерном пространстве. Решение линейных неравенств, определяющих внутреннюю область треугольника. Определение пределов, производных; исследование функции.
контрольная работа [567,1 K], добавлен 21.05.2013Определение дифференциальных уравнений в частных производных параболического типа. Приведение уравнения второго порядка к каноническому виду. Принцип построения разностных схем. Конечно-разностный метод решения задач. Двусторонний метод аппроксимации.
дипломная работа [603,8 K], добавлен 24.01.2013Матричный метод решения систем линейных алгебраических уравнений с ненулевым определителем. Примеры вычисления определителя матрицы. Блок-схема программы, описание объектов. Графический интерфейс, представляющий собой стандартный набор компонентов Delphi.
курсовая работа [1,4 M], добавлен 29.06.2014Обобщенные решения линейных дифференциальных уравнений. Основные примеры построения фундаментальных решений линейных дифференциальных операторов с постоянными коэффициентами, метод преобразования Фурье. Преимущества использования методов спуска.
курсовая работа [1,1 M], добавлен 10.04.2014Математическое объяснение метода Эйлера, исправленный и модифицированный методы. Блок-схемы алгоритмов, описание, текст и результаты работы программы. Решение обыкновенных дифференциальных (нелинейных) уравнений первого порядка с начальными данными.
курсовая работа [78,1 K], добавлен 12.06.2010
