Исследование метода Гаусса для плотных матриц с выбором главного элемента по столбцу

Министерство  образования и науки Российской Федерации

НОВОСИБИРСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ 

Кафедра экономической информатики 
 
 
 
 
 
 
 
 
 
 
 
 

Курсовой  проект

По курсу: «Численные методы»

Тема: «Исследование  метода Гаусса для плотных матриц с выбором главного элемента по столбцу» 
 
 
 
 
 
 
 
 
 

Выполнила студентка ОТЗ-707

Чучуйко Е.В. 

Научный руководитель:

К.т.н., доцент Сарычева О.М. 
 
 
 
 
 
 
 
 
 

Новосибирск  2009

 

      Содержание 

     1 Постановка задачи (математическое описание метода)

     2 Описание ПО

     3 Описание тестовых задач

     4 Анализ результатов счета

     5 Выводы

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

     7 ПРИЛОЖЕНИЕ 

 

      Постановка задачи (математическое описание метода) 

     Одним из самых распространенных методов  решения систем линейных уравнений  является метод Гаусса.Этот метод (который также называют методом последовательного исключения неизвестных) известен в различных вариантах уже более 2000 лет. 

     Вычисления  с помощью метода Гаусса заключаются  в последовательном исключении неизвестных  из системы для преобразования ее к эквивалентной системе с верхней треугольной матрицей. Вычисления значений неизвестных производят на этапе обратного хода. 

       Схема единственного деления.  Рассмотрим сначала простейший  вариант метода Гаусса, называемый  схемой единственного деления. 

     Прямой  ход состоит из n- 1 шагов исключения. 

     1-й  шаг. Целью этого шага является  исключение неизвестного x1 из уравнений  с номерами i = 2, 3, …, n. Предположим,  что коэффициент a11¹ 0. Будем  называть его главным элементом  1-го шага. 

     Найдем  величины 

     qi1 = ai1/a11  (i = 2, 3, …, n), 

     называемые  множителями 1-го шага. Вычтем последовательно  из второго, третьего, …, n-го уравнений  системы первое уравнение, умноженное соответственно на q21, q31, …, qn1. Это позволит обратить в нуль коэффициенты при x1 во всех уравнениях, кроме первого. В результате получим эквивалентную систему 

     a11x1 + a12x2 + a13x3 + … + a1nxn= b1 , 

     a22(1)x2 + a23(1)x3 + … + a2n(1)xn= b2(1) , 

     a32(1)x2 + a33(1)x3 + … + a3n(1)xn= b3(1) , 

     .  .  .  .  .  .  .  .  .  .  .  .  .  .  . 

     an2(1)x2 + an3(1)x3 + … + ann(1)xn= bn(1) . 

     в которой aij(1)и bij(1) вычисляются по формулам 

     aij(1) = aij − qi1a1j   ,   bi(1) = bi −  qi1b1. 

     2-й  шаг. Целью этого шага является  ислючение неизвестного x2 из уравнений  с номерами i = 3, 4, …, n. Пусть a22(1) ≠  0, где a22(1) – коэффициент, называемый главным (или ведущим) элементом 2-го шага. Вычислим множители 2-го шага 

     qi2 = ai2(1) / a22(1)  (i = 3, 4, …, n) 

     и вычтем последовательно из третьего, четвертого, …, n-го уравнения системы  второе уравнение, умноженное соответственно на q32, q42, …, qm2. В результате получим систему 

            a11x1+a12x2+   a13x3+  … + a1nxn =    b1   , 

                  a22(1)x2 +   a23(1)x3 +      …  +      a2n(1) =    b2(1)     , 

                          a33(2)x3 +… + a3n(2)xn =   b3(2)  , 

     .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  . 

                          an3(2)x3 +      … +      ann(2)xn= bn(2)     . 

     Здесь коэффициенты aij(2)и bij(2) вычисляются  по формулам 

     aij(2) = aij(1) – qi2a2j(1)  ,   bi(2) = bi(1) –  qi2b2(1). 

     Аналогично  проводятся остальные шаги. Опишем очередной k-й шаг. 

     k-й  шаг. В предположении, что главный  (ведущий) элемент k-го шага akk(k–1) отличен от нуля, вычислим множители  k-гошага 

     qik = aik(k–1) / akk(k–1)  (i = k + 1, …, n) 

     и вычтем последовательно из (k + 1)-го, …, n-го уравнений полученной на предыдущем шаге системы k-e уравнение, умноженное соответственно на qk+1,k, qk+2,k, …, qnk. 

     После (n - 1)-го шага исключения получим систему  уравнений 

            a11x1 +a12x2 +    a13x3 +  … +  a1nxn =     b1     , 

                  a22(1)x2 +  a23(1)x3 +      … +       a2n(1)xn =       b2(1)     , 

                          a33(2)x3 +      … +       a3n(2)xn =       b3(2)     , 

            .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  .  . 

                                      ann(n–1)xn =  bn(n–1). 

     матрица A(n-1) которой является верхней треугольной. На этом вычисления прямого хода заканчиваются. 

     Обратный  ход. Из последнего уравнения системы  находим xn. Подставляя найденное значение xn в предпоследнее уравнение, получим xn–1. Осуществляя обратную подстановку, далее последовательно находим xn–1, xn–2, …, x1. Вычисления неизвестных здесь проводятся по формулам 

     xn = bn(n–1) / ann(n–1), 

     xk = (bn(k–1) – ak,k+1(k–1)xk+1 – … – akn(k–1)xn) / akk(k–1), (k = n – 1, …, 1). 

     Необходимость выбора главных элементов. Заметим, что вычисление множителей, а также обратная подстановка требуют деления на главные элементы akk(k–1). Поэтому если один из главных элементов оказывыется равным нулю, то схема единственного деления не может быть реализована. Здравый смысл подсказывает, что и в ситуации, когда все главные элементы отличны от нуля, но среди них есть близкие к нулю, возможен неконтролируемый рост погрешности. 

     1.1.2. Метод Гаусса с выбором главного  элемента по столбцу (схема  частичного выбора). Описание метода. На k-м шаге прямого хода коэффициенты уравнений системы с номерами i = k + 1, …, nпреобразуются по формулам 

     aij(k) = aij(k–1)− qikakj , bi(k) = bi(k–1) − qikbk(k–1), i = k + 1, …, n. 

     Интуитивно  ясно, что во избежание сильного роста коэффициентов системы и связанных с этим ошибок нельзя допускать появления больших множителей qik. 

     В методе Гаусса с выбором главного элементоа по столбцу гарантируется, что |qik| ≤ 1 для всех k = 1, 2, …, n – 1 и i = k + 1, …, n. Отличие этого варианта метода Гаусса от схемы единственного деления заключается в том, что на k-м шаге исключения в качестве главного элемента выбирают максимальный по модулю коэффициент aikkпри неизвестной xkв уравнениях с номерами i = k + 1, …, n. Затем соответствующее выбранному коэффициенту уравнение с номером ik меняют местами с k-м уравнением системы для того, чтобы главный элемент занял место коэффициента 

akk(k-1). После этой перестановки исключение  неизвестного xk производят, как в схеме единственного деления. 

     Описание  ПО 

     Основной  особенностью языка MATLAB является его  широкие возможности по работе с  матрицами. Язык MATLAB является высокоуровневым  интерпретируемым языком программирования, включающим основанные на матрицах структуры  данных, широкий спектр функций, интегрированную среду разработки, объектно-ориентированные возможности и интерфейсы к программам, написанным на других языках программирования. Для работы использовалась среда MatLab 7.0 

     Описание  тестовых задач

В качестве проверки правильности работы алгоритма выбрана система уравнений

5х₁+х₂+х₃=5

Х₁+4х₂+х₃=4

Х₁+х₂+3х₃=3

Решения системы, полученные в ручную и в  программном режиме совпадают и  равны х₁= 0,76; х₂= 0,68; х₃=0,52_; 

     Анализ  результатов счета

       

     

       
 

     

     Выводы

     При увеличении размерности матрицы появляются ошибки вычислений. 

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

1. Калиткин Н.Н. Численные методы. - М.: Наука, 1978.
2. Н.В.Копченова, И.А. Марон. Вычислительная математика в примерах и задачах. - М.: Наука, 1972.
3. О.М. Сарычева Численные методы в экономике, конспект лекций, Новосибирск, 1995

 

 

      ПРИЛОЖЕНИЕ 

     main.,m 

     clc;

     r=10;

      

     %Исследование  размерности и обусловленности

     er1=[]; er2=[]; nn=[];

     eob1=[]; eob2=[];

      

     for n=2:r

         [A,b]=matrix(r);

         As=A;

         bs=b;

         A=[A b];

         A=gauss(A);

         [A, b]=AandB(A);

         x=obrxod(A, b);

         e=epsilon(As, bs, x);

         [E1, E2]=oshibka(e);

         er1=[er1 E1];

         er2=[er2 E2];

        

     %Исследование  обусловленности

         A=obusl(As);

         Apo=A;

         A=[A bs];

         [A]=gauss(A);

         [A, b]=AandB(A);

         [x]=obrxod(A, b);

         [e]=epsilon(Apo, bs, x);

         [E11, E22]=oshibka(e);

         eob1=[eob1 E11];

         eob2=[eob2 E22];

         nn=[nn n];

     end

      

     plot(nn, er1);

     title('Srednekvadratichnaya oshibka po razmernosti');

     pause;

      

     plot(nn, er2);

     title('Maksimalnaya oshibka po razmernosti');

     pause;

      

     plot(nn, eob1);

     title('Srednekvadratichnaya oshibka dlya plohoy obuslovlennosti');