Численные методы решения систем плохо обусловленных уравнений

Уже при k>3 метод сходится совершенно к другому решению. Зато на матрице Гильберта метод Гаусса-Зейделя дает очень хорошие результаты. Сравнивая его с методом Якоби на этой матрице (табл.

8 и табл.

6) видно, что метод Гаусса-Зейделя оказался менее чувствительным к спектральному радиусу матрицы СЛАУ.

6)Метод сопряженных градиентов.

Так как обе рассматриваемые матрицы являются симметричными, то из проекционных методов Крыловского типа был выбран MSG. Приведем результаты тестирования.

k (x*-x)/—x*— —b-A*x—/—b— iter 1 3.11097e-015 5.12002e-015 5 2 1.12024e-013 7.76077e-014 5 3 1.93552e-013 3.62478e-013 5 4 1.7421e-011 3.31552e-012 5 5 1.34447e-010 9.73895e-011 5 6 3.03046e-010 3.79453e-010 5 7 0.603 023 2.81439e-009 4 8 0.603 023 2.81439e-010 4 9 0.603 023 2.8144e-011 4 10 0.603 023 2.81444e-012 4 11 0.603 023 2.81514e-013 4 12 0.603 023 2.81967e-014 4 13 0.603 023 2.7274e-015 4 14 0.603 023 2.30459e-016 4 15 0.603 023 1.67684e-016 4 16 0.603 023 2.18275e-016 4 17 0.603 023 2.18275e-016 4 18 0.603 023 2.18275e-016 4 19 0.603 023 2.18275e-016 4 Табл.

9 Метод сопряженных градиентов (матрица первого типа)

k (x*-x)/—x*— —b-A*x—/—b— iter 1 0,00E+00 0,00E+00 0 2 4.96507e-016 4.79494e-016 2 3 1.45408e-013 1.98275e-013 3 4 1.50616e-009 2.23171e-009 4 5 3.41715e-008 7.02783e-009 6 6 0.305 777 1.74313e-009 6 7 0.595 865 2.92454e-009 6 8 0.140 773 2.91885e-010 8 9 0.228 505 1.4111e-010 8 10 0.326 601 4.09149e-010 8 11 0.432 177 1.04466e-009 8 12 0.542 998 2.25748e-009 8 13 0.657 355 4.32574e-009 8 14 0.773 947 7.58103e-009 8 15 0.223 795 1.29039e-009 10 16 0.264 998 1.94236e-009 9 17 0.3 077 7.73266e-010 10 18 0.351 611 1.33605e-009 10 19 0.396 485 8.76969e-010 10 Табл.

10 Метод сопряженных градиентов (матрица Гильберта)

Из табл.

9 видно, что MSG для k < 7 дает очень точный результат, но при большем k, метод начинает сходиться не к истинному решению. Это связанно с тем, что с возрастом k начинает нарушаться положительная определенность матрицы A, что и влечет за собой ошибку в решении. Но в случае матрицы Гильберта (табл.

10) MSG показывает достаточно хороший результат. Заметим, что число итераций в MSG, как и должно быть, не превышает размерности матрицы СЛАУ, что является показателем того, что метод действительно нашел решение (только не всегда истинное, как в случае исследования на матрицах первого типа).

7) Сравнение алгоритмов.

Сравнивая результаты исследования на матрицах первого типа (табл.

1, табл.

3, табл.

5, табл.

7, табл.

9) видно, что наиболее лучший результат показали прямые методы, а хуже всех себя зарекомендовали классические итерационные методы. MCG на этой матрице хоть и сходился к решению за очень маленькое число итераций, при k>6 сходился уже не к точному решению. Тестирование на таблице Гильберта показали совершенно другие результаты (табл.

2, табл.

4, табл.

6, табл.

8, табл.

10). Самым устойчивым в этом случае оказался MSG. Наиболее не приспособленным опять же оказался метод Якоби. Не смотря на то, что при размерности матрицы Гильберта k < 10 прямые методы показали самые хорошие результаты, при k>10 эти методы давали очень большую погрешность.

Заключение

В дано работе приведены и реализованы основные методы решения СЛАУ. Исследование влияния числа обусловленности матрицы на результаты методов показали, что в различных случаях могут применяться различные методы и нельзя подобрать универсальный метод решения СЛАУ. В каждом конкретном случае (в зависимости от задачи, от свойств матрицы) нужно подбирать решатель способный достаточно хорошо справится с поставленной задачей. В работе также были представлены методы предобусловливания, которые также позволяют бороться с числом обусловленности. Помимо диагонального предобусловливания и LU существуют и более сложные методы, которые строятся уже исключительно из самой задачи, а не ее матрицы СЛАУ. К таким методам относятся различные многосеточные методы, используемые в решении дифференциальных уравнений второго порядка методом конечных элементов или конечных объемов[6].

Список литературы

1 Ильин В. П. Методы неполной факторизации для решения линейных систем.

М.:Физматлит, 1995.

2 Баландин М. Ю., Шурина Э. П. Методы решения СЛАУ большой размерности: Учеб. Пособие. — Новосибирск: Изд-во НГТУ, 2000.

3 Трофимов В. П. Численные методы решения систем линейных алгебраических уравнений.

Воронеж, 2004.

4 Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы.

М.: Лаборатория Базовых Знаний, 2002.

5 Курцева К. П., Шешуков Е. Г. Методы вычислений. Часть 1. Численные методы алгебры.

Зеленодольск, 2007.

6 Рояк М. Э., Соловейчик Ю. Г., Шурина Э. П. Сеточные методы решения кравевых задач математической физики: Учеб. Пособие. — Новосибирск: Изд-во НГТУ, 1998.

7 Калиткин Н. Н. Численные методы.

М.: Наука, 1978.

Приложение 1

Выбрать и начальное приближение

Создать матрицу Создать вектор Начало

Для

Для

Увеличить

Увеличить

Решить СЛАУ

Конец Приложение 1. Алгоритм GMRES

Приложение 2

Текст программы:

#pragma once

#include

#include

#include

#include

#include

#include

#include

using namespace std;

//Метод гаусса с выбором ведущего элемента

void GAUSS (double**A, double*x, double*b, int dims)

{

int i, j, k;

double max;

int nmax;

//прямой ход

for (i = 0; i < dims — 1; i++)

{

max = fabs (A[i][i]);

nmax = i;

for (j = i + 1; j < dims; j++)

if (max < fabs (A[j][i]))

{

nmax = j;

max = fabs (A[j][i]);

}

if (nmax ≠ i)

{

for (j = i; j < dims; j++)

{

max = A[i][j];

A[i][j] = A[nmax][j];

A[nmax][j] = max;

}

max = b[i];

b[i] = b[nmax];

b[nmax] = max;

}

for (j = i + 1; j < dims; j++)

{

for (k = i + 1; k < dims; k++)

A[j][k] -= A[i][k] / A[i][i] * A[j][i];

b[j] -= b[i] / A[i][i] * A[j][i];

A[j][i] = 0.0;

}

}

//Обратный ход

for (i = dims — 1; i >= 0; i—)

{

for (j = dims — 1; j > i; j—)

b[i] -= A[i][j] * x[j];

x[i] = b[i] / A[i][i];

}

};

//умножение матрицы на вектор

void Mult (double **A, double *x, double*b, int dims)

{

int i, j;

for (i = 0; i < dims; i++)

{

b[i] = 0.0;

for (j = 0; j < dims; j++)

b[i] += A[i][j] * x[j];

}

};

//Генератор матрицы (1)

void Gen₁(double**A, int k, int dims)

{

A[1][0] = -0.0;

A[2][0] = -3.0;

A[2][1] = -1.0;

A[3][0] = -4.0;

A[3][1] = -2.0;

A[3][2] = -3.0;

A[4][0] = -0.0;

A[4][1] = -2.0;

A[4][2] = -4.0;

A[4][3] = -1.0;

int i, j;

double sum;

for (i = 0; i < dims; i++)

{

for (j = i + 1; j < dims; j++)

{

A[i][j] = A[j][i];

}

}

for (i = 0; i < dims; i++)

{

sum = 0;

for (j = 0; j < dims; j++)

if (j ≠ i)

{

sum += A[i][j];

}

A[i][i] = - sum;

}

sum = 1;

for (i = 0; i < k; i++)

{

sum = sum / 10;

}

A[0][0] += sum;

};

//подсчет нормы

double norma (double *x, int dims)

{

int i;

double temp = 0.0;

for (i = 0; i < dims; i++)

{

temp += x[i] * x[i];

}

return sqrt (temp);

};

//Генератор матрицы Гильберта

void Gen₂(double**A, int k, int dims)

{

int i, j;

for (i = 0; i < k; i++)

for (j = 0; j < k; j++)

A[i][j] = 1.0 / (i + j + 1.0);

};

//Метод Якоби

int Jakobi (double**A, double *x, double *b, int dims)

{

int i, j, k;

double *x0;

x0 = new double[dims];

for (i = 0; i < dims; i++)

x0[i] = 0.0;

double *nevjazka;

nevjazka = new double[dims];

for (i = 0; i < dims; i++)

{

nevjazka[i] = b[i];

}

double norm_b = norma (b, dims);

int max_iter = 200;

double gamma = 1.0e-7;

double sum;

for (i = 0; i < max_iter && norma (nevjazka, dims) / norm_b >= gamma; i++)

{

for (j = 0; j < dims; j++)

{

sum = 0;

for (k = 0; k < dims; k++)

{

if (k ≠ j)

{

sum += A[j][k] * x0[k];

}

}

x[j] = 1.0 / A[j][j] * (b[j] - sum);

}

for (j = 0; j < dims; j++)

x0[j] = x[j];

Mult (A, x0, nevjazka, dims);

for (j = 0; j < dims; j++)

nevjazka[j] = b[j] - nevjazka[j];

cout << «iter = «<<» nev = «<< norma (nevjazka, dims) / norm_b << endl;

}

return i;

};

//Метод Гаусса-Зейделя

int Gaus_Zeidel (double**A, double*x, double*b, int dims)

{

int i, j, k;

double *x0;

x0 = new double[dims];

for (i = 0; i < dims; i++)

x0[i] = 0.0;

double *nevjazka;

nevjazka = new double[dims];

for (i = 0; i < dims; i++)

{

nevjazka[i] = b[i];

}

double norm_b = norma (b, dims);

int max_iter = 20 000;

double gamma = 1.0e-7;

double sum;

for (i = 0; i < max_iter && norma (nevjazka, dims) / norm_b >= gamma; i++)

{

for (j = 0; j < dims; j++)

{

sum = 0;

for (k = 0; k < j; k++)

{

sum += A[j][k] * x[k];

}

for (k = j + 1; k < dims; k++)

{

sum += A[j][k] * x0[k];

}

x[j] = 1.0 / A[j][j] * (b[j] - sum);

}

for (j = 0; j < dims; j++)

x0[j] = x[j];

Mult (A, x0, nevjazka, dims);

for (j = 0; j < dims; j++)

nevjazka[j] = b[j] - nevjazka[j];

cout << «iter = «<<» nev = «<< norma (nevjazka, dims) / norm_b << endl;

}

return i;

};

//Метод LU-разложения

void LU (double **A, double*x, double*b, int dims)

{

int i, j, k;

double sum;

for (i = 0; i < dims; i++)

{

for (j = 0; j < i; j++)

{

sum = 0.0;

for (k = 0; k < j; k++)

{

sum += A[i][k] * A[k][j];

}

A[i][j] = (A[i][j] - sum) / A[j][j];

}

for (j = 0; j <= i; j++)

{

sum = 0.0;

for (k = 0; k < j; k++)

{

sum += A[k][i] * A[j][k];

}

A[j][i] = A[j][i] - sum;

}

}

double *y;

y = new double[dims];

//L^(-1)y=b

for (i = 0; i < dims; i++)

{

sum = 0;

for (j = 0; j < i; j++)

sum += A[i][j] * y[j];

y[i] = b[i] - sum;

}

//U^(-1)x=y

for (i = dims — 1; i >= 0; i—)

{

sum = 0.0;

for (j = i + 1; j < dims; j++)

sum += A[i][j] * x[j];

x[i] = (y[i] - sum) / A[i][i];

}

};

//Скалярное умножение векторов

double Scalar (double *v, double *u, int dims)

{

double temp = 0;

int i;

for (i = 0; i < dims; i++)

temp += v[i] * u[i];

return temp;

};

//Метод MSG

int MSG (double**A, double *x, double*b, int dims)

{

double alfa, betta;

double *r0, *r;

double *p, *temp;

int i, max_iter, j;

double gamma = 1.0e-8;

max_iter = 10 000;

r0 = new double [dims];

r = new double [dims];

p = new double [dims];

temp = new double [dims];

for (i = 0; i < dims; i++)

x[i] = 1.0;

Mult (A, x, r0, dims);

for (i = 0; i < dims; i++)

{

p[i] = r0[i] = b[i] - r0[i];

}

for (i = 0; i < max_iter && norma (r0, dims) / norma (b, dims) > gamma; i++)

{

Mult (A, p, temp, dims);

alfa = Scalar (r0, r0, dims) / Scalar (temp, p, dims);

for (j = 0; j < dims; j++)

{

x[j] = x[j] + alfa * p[j];

}

for (j = 0; j < dims; j++)

{

r[j] = r0[j] - alfa * temp[j];

}

betta = Scalar (r, r, dims) / Scalar (r0, r0, dims);

for (j = 0; j < dims; j++)

{

p[j] = r[j] + betta * p[j];

r0[j] = r[j];

}

cout << norma (r, dims) / norma (b, dims) << endl;

}

return i;

};

//Главная часть

void main ()

{

double norma_b;

double **A;

double **A2;

double *b;

double *b2;

double *x;

double *rezult;

double *pogr;

double *nevjazka;

int dims = 5;

int i, j, k;

ofstream f («output₂_MSG.txt»);

for (k = 1; k < 20; k++)

{

dims = k;

A = new double* [dims];

for (i = 0; i < dims; i++)

A[i] = new double[dims];

A2 = new double* [dims];

for (i = 0; i < dims; i++)

A2[i] = new double[dims];

b = new double[dims];

b2 = new double[dims];

x = new double[dims];

rezult = new double[dims];

pogr = new double[dims];

nevjazka = new double[dims];

for (i = 0; i < dims; i++)

{

x[i] = i + 1;

// rezult[i] = 1;

}

Gen₂(A, k, dims);

Gen₂(A2, k, dims);

Mult (A, x, b, dims);

Mult (A2, x, b2, dims);

norma_b = norma (b2, dims);

// GAUSS (A, rezult, b, dims);

// LU (A, rezult, b, dims);

// j = Gaus_Zeidel (A, rezult, b, dims);

// j = Jakobi (A, rezult, b, dims);

j = MSG (A, rezult, b, dims);

f << k << endl;

f <<" iter = «<

for (i = 0; i < dims; i++)

pogr[i] = rezult[i] - x[i];

f <<" pogr_otn = «<<

Mult (A2, rezult, nevjazka, dims);

for (i = 0; i < dims; i++)

nevjazka[i] = b2[i] - nevjazka[i];

f <<" nevjazka_otn = «<<

f << endl;

delete A;

delete A2;

delete b;

delete b2;

delete x;

delete rezult;

delete pogr;

delete nevjazka;

}

f.close ();

// getchar ();

};