КИЇВСЬКИЙ НАЦІОНАЛЬНИЙ УНІВЕРСИТЕТ ІМЕНІ ТАРАСА ШЕВЧЕНКА Лабораторна робота 2
Зміст: Постановка задачі. Обґрунтування методів. Результати. Код програми. Висновок. 1. Постановка задачі. Знайти розв’язок системи лінійних рівнянь Ax=b, де А: b: |1 2 0 0 0 0 0 0 …| |1| |1 3 2 0 0 0 0 0 …| |2| |0 1 3 2 0 0 0 0 …| |3| |0 0 1 3 2 0 0 0 …| |4| |0 0 0 1 3 2 0 0 …| |5| |0 0 0 0 1 3 2 0 …| |6| |0 0 0 0 0 1 3 2 …| |7| |0 0 0 0 0 0 1 3 …| |8| |…………………| |..| Методами Якобі та Прогонки Cx=b C: b: |1 1 0 0 0 0 0 0 …| |1| |a1 b1 1 0 0 0 0 0 …| |2| |0 a2 b2 1 0 0 0 0 …| |3| |0 0 a3 b3 1 0 0 0 …| |4| |0 0 0 a4 b4 1 0 0 …| |5| |0 0 0 0 a5 b5 1 0 …| |6| |0 0 0 0 0 a6 b6 1 …| |7| |0 0 0 0 0 0 a7 b7 …| |8| |…………………………...| |..| (аі = 2+2/х; bi = 1/x) Методами Зейделя та Прогонки Для кожної системи знайти нев'язку та Розв’язати з точністю до 0,001 і визначити кількість ітерацій. 2. Обґрунтування методів. Метод прогонки Метод прогонки или алгоритм Томаса (англ. Thomas algorithm) используется для решения систем линейных уравнений вида /, где A — трёхдиагональная матрица. Описание метода Система уравнений / равносильна соотношению / Метод прогонки основывается на предположении, что искомые неизвестные связаны рекуррентным соотношением: / где / Используя это соотношение, выразим xi-1 и xi через xi+1 и подставим в уравнение (1): /, где Fi — правая часть i-го уравнения. Это соотношение будет выполняться независимо от решения, если потребовать / Отсюда следует: / Из первого уравнения получим: / После нахождения прогоночных коэффициентов ? и ?, используя уравнение (2), получим решение системы. При этом, / / / Другим способом объяснения существа метода прогонки, более близким к терминологии конечно-разностных методов и объясняющим происхождение его названия, является следующий: преобразуем уравнение (1) к эквивалентному ему уравнению / c наддиагональной матрицей /. Вычисления проводятся в два этапа. На первом этапе вычисляются компоненты матрицы / и вектора /, начиная с / до / / / и / / На втором этапе, для / вычисляется решение: / / Такая схема вычисления объясняет также английский термин этого метода «shuttle». Для применимости формул метода прогонки достаточно свойства диагонального преобладания у матрицы A.
Метод Якоби Возьмём систему линейных уравнений: /, где / Или / Описание метода Для того, чтобы построить итеративную процедуру метода Якоби, необходимо провести предварительное преобразование системы уравнений / к итерационному виду /. Оно может быть осуществлено по одному из следующих правил: / / / где в принятых обозначениях D означает матрицу, у которой на главной диагонали стоят соответствующие элементы матрицы A, а все остальные нули; тогда как матрицы U и L содержат верхнюю и нижнюю треугольные части A, на главной диагонали которых нули, E — единичная матрица. Тогда процедура нахождения решения имеет вид: / где k счётчик итерации. В отличие от метода Гаусса-Зейделя мы не можем заменять / на / в процессе итерационной процедуры, т.к. эти значения понадобятся для остальных вычислений. Это наиболее значимое различие между методом Якоби и методом Гаусса-Зейделя решения СЛАУ. Таким образом на каждой итерации придётся хранить оба вектора приближений: старый и новый. Условие сходимости Приведем достаточное условие сходимости метода. / Теорема.Пусть /. Тогда при любом выборе начального приближения /: метод сходится; скорость сходимости метода равна скорости сходимости геометрической прогрессии со знаменателем /; верна оценка погрешности: /.
Условие остановки Условие окончания итерационного процесса при достижении точности / в упрощённой форме имеет вид: / (Существует более точное условие окончания итерационного процесса, которое более сложно и требует дополнительных вычислений)
Метод Зейделя Возьмём систему: /, где / Или / И покажем, как её можно решить с использованием метода Гаусса-Зейделя. Метод Чтобы пояснить суть метода, перепишем задачу в виде: / Здесь в j-м уравнении мы перенесли в правую часть все члены, содержащие xi , для i > j. Эта запись может быть представлена: / где в принятых обозначениях D означает матрицу, у которой на главной диагонали стоят соответствующие элементы матрицы A, а все остальные нули; тогда как матрицы U и L содержат верхнюю и нижнюю треугольные части A, на главной диагонали которых нули. Итерационный процесс в методе Гаусса-Зейделя строится по формуле / после выбора соответствующего начального приближения /. Метод Гаусса-Зейделя можно рассматривать как модификацию метода Якоби. Основная идея модификации состоит в том, что новые значения / используются здесь сразу же по мере получения, в то время как в методе Якоби они не используются до следующей итерации: / где / Таким образом i-тая компонента (k + 1)-го приближения вычисляется по формуле: / Условие сходимости Приведём достаточное условие сходимости метода. / Теорема.Пусть /, где / – матрица, обратная к /. Тогда при любом выборе начального приближения /: метод Гаусса-Зейделя сходится; скорость сходимости метода равна скорости сходимости геометрической прогрессии со знаменателем /; верна оценка погрешности: /.
Условие окончания Условие окончания итерационного процесса Зейделя при достижении точности / в упрощённой форме имеет вид: / Более точное условие окончания итерационного процесса имеет вид / и требует больше вычислений. Хорошо подходит для разреженных матриц.
Підрахування констант:
3. Результати. 4. Код програми. #include <iostream> #include <math.h> #define MaxN 100 #define eps 0.001 using namespace std; /// our matrix double A[MaxN][MaxN]; /// vector x - we must find it double x[MaxN]; /// first answer double x0[MaxN]; /// vector of free elements double b[MaxN]; /// nevyazka double nev[MaxN]; /// matrix range int n = 0; ///have an answer bool haveanswer = false; ///second - ??? double second = 0; /// init matrix for first lab void Init1() { for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { A[i][j] = 0; } A[i][i] = 3; if (i == 0) { A[0][0] = 1; } else { A[i][i-1] = 1; } if (i <= n-1) { A[i][i+1] = 2; } b[i] = i + 1; } haveanswer = false; } ///help functions double a(double x) { return 2+2/x; } double p(double x) { return 1/x; } /// init matrix for second lab void Init2() { for (int i = 0; i < n; i++) { for (int j = 0; j < n; j++) { A[i][j] = 0; } A[i][i] = a(i+1); if (i > 0) { A[i][i-1] = p(i); } if (i <= n-1) { A[i][i+1] = 1; } b[i] = i + 1; } haveanswer = false; } ///print Ax=b; /// if we don't know x than print xi else numbers void print() { for (int i = 0; i < n; i++) { cout << "| "; for (int j = 0; j < n; j++) { cout << A[i][j] << " "; } cout << (i==(int)(n/2)?"|x|":"| |") << (haveanswer?"":"x") << (haveanswer?x[i]:i) << (i==(int)(n/2)?"|=|":"| |") << b[i] << "|" << endl; } } /// reset vector x void ResetX() { for(int i = 0; i < n; i++) { x[i] = 0; } } /// method Jacobi void Jacobi() { cout << "Method Jacobi starts..." << endl; double norm; do { for (int i = 0; i < n; i++) { x0[i] =- b[i]; for (int g = 0; g < n; g++) { if (i != g) x0[i] += A[i][g] * x[g]; } x0[i] /= -A[i][i]; } norm = fabs(x[0] - x0[0]); for (int h = 0; h < n; h++) { if (fabs(x[h] - x0[h]) > norm) norm = fabs(x[h] - x0[h]); x[h] = x0[h]; } } while (norm > eps); haveanswer = true; } /// method Sweep void Sweep() { cout << "Method Sweep starts..." << endl; /// 3-diagonal to 2-diagonal A[0][1] = A[0][1] / A[0][0]; b[0] = b[0] / A[0][0]; A[0][0] = 1; for(int i = 1; i < n-1; i++) { A[i][i+1] = A[i][i+1] / (A[i][i] - A[i][i-1] * A[i-1][i]); b[i] = (b[i] - A[i][i-1] * b[i-1])/(A[i][i] - A[i][i-1] * A[i-1][i]); A[i][i-1] = 0; A[i][i] = 1; } b[n-1] = (b[n-1] - A[n-1][n-2] * b[n-2])/(A[n-1][n-1] - A[n-1][n-2] * A[n-2][n-1]); A[n-1][n-2] = 0; A[n-1][n-1] = 1; /// find answer: x[n-1] = b[n-1]; for(int i = n-2; i >= 0; i--) { x[i] = b[i] - x[i+1] * A[i][i+1]; } haveanswer = true; } ///help function bool converge(double *xk, double* xkp) { for (int i = 0; i < n; i++) { if (fabs(xk[i]-xkp[i]) >= eps) return false; } return true; } /// method Zeidel void Zeidel() { cout << "Method Zeidel starts..." << endl; do { for(int i = 0; i < n; i++) { double var = 0; for(int j = 0; j < n; j++) if(j != i) var += (A[i][j]*x[j]); x0[i] = x[i]; x[i]=(b[i] - var)/A[i][i]; } } while(!converge(x,x0)); haveanswer = true; } /// print nevyazka void Print_Nevyazka() { cout << "Nevyazka:" << endl; for(int i = 0; i < n; i++) { nev[i] = 0; for(int j = 0; j < n; j++) { nev[i] += A[i][j] * x[i]; } nev[i] = b[i] - nev[i]; cout << "| " << nev[i] << " |" << endl; } } /// print don't remember what :( void Print_Second() { cout << "Second:" << endl; /// calculate it /// print it cout << second << end; } /// main function int main() { cout << "Input n:" << endl; cin >> n; Init1(); print(); ResetX(); Jacobi(); print(); Init1(); Print_Nevyazka(); Print_Second(); ResetX(); Sweep(); print(); Init1(); Print_Nevyazka(); Print_Second(); cout<<"-----------------------------------------------------"<<endl; Init2(); print(); ResetX(); Zeidel(); print(); Print_Nevyazka(); Print_Second(); Init2(); ResetX(); Sweep(); print(); Init2(); Print_Nevyazka(); Print_Second(); system("pause"); return 0; } 5. Висновок.
Методи зійшлися, отже зроблені вірно. Похибка для таких методів апріорі відома: 0.001 – адже заздалегідь ми задаємо точність підрахунків.