Ирина Эланс

Автор который поможет с любыми образовательными и учебными заданиями

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

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

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

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

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

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

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

Выполнила студентка ОТЗ-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_;

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

Выводы

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

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

Калиткин Н.Н. Численные методы. - М.: Наука, 1978.
Н.В.Копченова, И.А. Марон. Вычислительная математика в примерах и задачах. - М.: Наука, 1972.
О.М. Сарычева Численные методы в экономике, конспект лекций, Новосибирск, 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');

pause;

plot(nn, eob2);

title('Maksimalnaya oshibka dlya plohoy obuslovlennosti');

pause;

%Исследование разряженной матрицы

A=[2 0 0 1;0 1 0 0; 0 0 1 0; 1 0 0 3];

b=[1; 3; 5; 2];

A=As;

A=[A bs];

[A]=gauss(A);

[A, b]=AandB(A);

[x]=obrxod(A, b);

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

[E1, E2]=oshibka(e)

gauss.m

function[A]=gauss(A)

n=size(A, 2);

for k=1:n-1

[A, k]=perest(A, k, n);

if k==n-1

break;

end

[A, k]=normir(A, k,n);

if k==n-1

break;

end

[A]=iskl(A, k, n);

end

perest.m

function[A, k]=perest(A, i, n)

ns=0;

j=1;

max=0;

while max==0

max=A(j, i);

for jn=j:n-1

if max<A(jn, i)

max=A(jn, i);

ns=jn;

end

if ns>0

for in=1:n

S=A(ns, in);

A(ns, in)=A(j, in);

A(j, in)=S;

end

k=i-1;

if i<n-1

k=i;

end

iskl.m

function[A]=iskl(A, i, n)

j=i+1;

jm=j-1;

for jn=j:n-1

d=-A(jn, i)/A(jm, i);

for in=i:n

A(jn, in)=A(jm, in)*d+A(jn, in);

end

normir.m

function[A, k]=normir(A, i, n)

j=i;

g=0;

while g==0

g=A(j, i);

if g==0

i=i+1;

end

if i==n-1

k=i;

else k=i;

end

for in=i:n

A(j, in)=A(j, in)/g;

end

obrxod.m

function[x]=obrxod(A, b)

s=0;

n=size(A, 1);

x=zeros(n, 1);

x(n)=b(n);

for is=1:n-1

i=n-is;

for ic=i+1:n

s=s+A(i, ic)*x(ic);

end

x(i)=b(i)-s;

s=0;

end

AandB.m

function[A, b]=AandB(Ad)

n=size(Ad, 1);

A=[];

for i=1:n

A=[A Ad(i, i)];

end

b=[Ad(:, n+1)];

epsilon.m

function[e]=epsilon(A, b, x)

n=size(A, 1);

e=zeros(n, 1);

for j=1:n

for i=1:n

e(j)=e(j)+A(j, i)*x;

end

e(j)=e(j)-b(j);

end

oshibka.m

function[E1, E2]=oshibka(e)

n=size(e, 1);

E1=0;

E2=zeros(1, n);

for i=1:n

E1=E1+e(i)^2;

end

E1=sqrt(E1);

for i=1:n

E2(i)=abs(e(i));

end

E2=max(E2);

obusl.m

function[A]=obusl(A)

A=A+10^2;

matrix.m

function[A, b]=matrix(n);

A=rand(n, n)*10;

b=rand(n, 1)*10;

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