Министерствообразования и науки Республики Беларусь
Белорусскийгосударственный университет
информатики ирадиоэлектроники
Факультетинформационных технологий и управления
Кафедра Вычислительных Методов и ПрограммированияПОЯСНИТЕЛЬНАЯ ЗАПИСКА
к курсовойработепо программированию
на тему:
«Решениесистемы линейных уравнений»
Выполнил: Принял:
ст.гр.020603 Навроцкий А.А.
Червоный А.В. Минск 2001г.
Содержание
Введение.
1. Анализсуществующих методов решения задачи.
2. Описаниеиспользуемого метода.
3. Анализрезультатов.
Вывод.
Список использованной литературы.
Приложение (распечатка программы, результатов).
ВведениеРешениесистем линейных алгебраических уравнений (СЛАУ) является одной из основныхзадач линейной алгебры. Эта задача имеет важное прикладное значение при решениинаучных и технических проблем. Кроме того, является вспомогательной приреализации многих алгоритмов вычислительной математики, математической физики,обработки результатов экспериментальных исследований.
Применяемые на практикечисленные методы решения СЛАУ делятся на две группы -прямые и итерационные.
В прямых (илиточных) методах решение системы получают за конечное число арифметическихдействий. К ним относятся известное правило Крамера нахождения решения спомощью определителей, метод последовательного исключения неизвестных (методГаусса) и его модификации, метод прогонки и другие. Сопоставление различныхпрямых методов проводится обычно по числу арифметический действий, необходимыхдля получения решения. Прямые методы являются универсальными и применяются длярешения систем до порядка 103. Отметим, что вследствие погрешностейокругления при решении задач на ЭВМ прямые методы на самом деле не приводят кточному решению системы.
Итерационные (или приближенные) методы являютсябесконечными и находят решение системы как предел при k®¥последовательных приближений x(k), где k — номер итерации. Обычно задается точность e, и вычисления проводятся до тех пор, пока не будет выполненаоценка ºx(k)– x(k-1)ºn(e), которое необходимо провести для получения заданнойточности, для многих методов можно найти из теоретических рассмотрений.Качество различных итерационных методов можно сравнивать по необходимому числуитераций n(e). Эти методы особенно предпочтительны для систем сматрицами специального вида — симметричными, трехдиагональными, ленточными ибольшими разреженными матрицами.
Выбор средыпрограммирования.
Послепроведенного обзора программных средств мы выбрали среду программированиянаиболее подходящую нам как очень удобное средство для разработки данногопрограммного продукта. DELPHI5.0 является наиболее выгодной нам средой программирования.
1. Анализ существующихметодов решения задачи
Прямые методы решенияСЛАУ. К прямым (илиточным) методам решения СЛАУ относятся алгоритмы, которые в предположении, чтовычисления ведутся без округлений, позволяют получить точное решение системы законечное число арифметических действий. Чаще всего решение задач такимиметодами осуществляется поэтапно: на первом этапе систему преобразуют к томуили иному простому виду, на втором — решают упрощенную систему и получаютзначения неизвестных.
Запишем систему линейныхалгебраических уравнений в развернутом виде:
/>
где x1, x2,...,xn — неизвестные величины, b1, b2,...,bn- элементы правой части. Если определитель системыотличен от нуля, то она имеет единственное решение. Для удобства дальнейшихпреобразований обозначим элементы правой части аi(n+1) и запишем расширенную матрицуразмерами n´(n+1), которая содержит всю информацию о системе:
A =/>.
С этой матрицей можнообращаться так же, как и с системой — переставлять строки, прибавлять кратноеодной строки к другой, исключая неизвестные и приводя матрицу к треугольномуили диагональному виду.
Приведем формальное описаниесхем некоторых прямых методов.
Метод Гаусса (схемаединственного деления). Алгоритм метода состоит из двух этапов. Первый этап называется прямымходом метода и заключается в последовательном исключении неизвестных изуравнений, т.е. в приведении матрицы А к верхнему треугольному виду(ниже главной диагонали все нули). Для этого на первом шаге разделим первоеуравнение системы на а11 (предположим, что коэффициента11¹0, в противном случае осуществляем перестановку уравнений системы). Обозначимкоэффициенты полученного приведенного уравнения />, домножим его на коэффициента21 и вычтем из второго уравнения системы, исключая тем самымх1 из второго уравнения (обнуляя коэффициента12матрицы). Поступим аналогично с остальными уравнениями и получим новую систему,матрица которой в первом столбце, кроме первого элемента, содержит только нули,т.е.
/>.
Первое уравнение вдальнейших преобразования не участвует. Описанный выше процесс исключениянеизвестных применим к матрице /> размерами(n-1) n. После k аналогичных шагов получим k приведенных уравнений с коэффициентами
/>
и матрицу /> размерами (n— k) (n— k+1), элементы которой вычисляются по формулам
/>.
Элементы />, на которые осуществляетсяделение, называются ведущими элементами метода Гаусса и не должныравняться нулю. Прямой ход метода Гаусса заканчивается после n шагов определением />.
Обратный ход метода Гаусса заключается впоследовательном определении компонент решения, начиная схnи заканчивая х1, по следующим формулам:
/>
Метод Гаусса с выборомглавного элемента. Метод заключается в том, что при прямом ходе в алгоритме метода Гаусса накаждом шаге исключения производится выбор наибольшего по модулю элемента вкачестве ведущего. Этого достигают перестановкой строк или столбцовматрицы коэффициентов. Наиболее распространённой в вычислительной практикеявляется стратегия выбора главного элемента столбца — нахождениемаксимального по модулю элемента k-гостолбца матрицы /> и использованиеего в качестве ведущего элемента на k-м шаге исключения. В этом случае для невырожденныхсистем гарантируется, что ведущие элементы не равны нулю, и уменьшается погрешностьпри делении и последующем вычитании при преобразованиях. Рекомендуется такжемасштабировать предварительно каждое уравнение исходной системы, разделив наего наибольший по абсолютной величине коэффициент. Это делает рост элементовпромежуточных матриц ограниченным.
Метод оптимальногоисключения.Вцелях экономии оперативной памяти (примерно в 4 раза) операции прямого иобратного хода метода Гаусса выполняются попеременно. На первом шаге послеприведения первого уравнения исключается неизвестное x1 извторого уравнения, а затем с помощью приведенного второго уравнения — неизвестноеx2 из первого. После (k-1) таких шагов матрица системы имеет вид
/>.
На k-м шаге, используя первыеk уравнений, исключаем неизвестные x1,..,xk из (k+1)-го уравнения. Затем посредством этого уравненияисключается неизвестное xk+1 из первых kуравнений и т.д. В результате прямогохода матрица системы приводится к диагональному виду с единицами на главнойдиагонали. При этом отпадает необходимость обратного хода, поскольку столбецправой части приведенной матрицы /> иявляется вектором решения.
Метод Гаусса-Жордана. Эта модификация метода Гауссанезначительно отличается от метода оптимального исключения. Операции исключенияпеременных для каждого приводимого уравнения осуществляют не только ниже, но ивыше главной диагонали. Операции с первым уравнением системы полностьюаналогичны стандартной схеме. Второе уравнение системы после приведения идомножения на соответствующие коэффициенты вычитаем не только из третьего ипоследующих уравнений, но и из первого. В результате k таких шагов получаем матрицу
/>.
Как и в методеоптимального исключения, матрица системы приводится к диагональному виду ивектором решения является столбец />.
LU -разложение. Матрицу A можно представить ввиде произведения нижней треугольной матрицы (включая диагональ) L (lower)и верхней треугольной матрицы U ( upper ), т.е. A=LU. Это равенство равносильно n2числовымравенствам
/>.
Разложение матрицы Aна множители обычно получают посредством алгоритма, который называется компактнойсхемой метода Гаусса. Элементыlim и Umi могут быть вычислены по формулам
/>
Тогда решение системы Ax=bсводится к последовательному решению двух систем — Ly=b и Ux=y.
Рассмотренный метод можноприменять к решению серии систем с одной и той же матрицей.
Метод простых итераций(Якоби).
Для решения итерационнымметодом система линейных алгебраических уравнений Ax= bдолжна быть приведена к виду x= Gx+f, где G — некоторая матрица, f-преобразованный вектор свободных членов. Затем выбирается начальное приближение- произвольный вектор x(0) — и строится рекуррентная последовательность векторов x(1), x(2),...,x(k),... по формуле
/>.
Для сходимости этойпоследовательности при любом начальном приближении необходимо и достаточно,чтобы все собственные значения матрицы G были по абсолютной величине меньше единицы. Напрактике это трудно проверить, и обычно пользуются достаточными условиямисходимости — итерации сходятся, если какая-нибудь норма матрицы меньшеединицы, т.е.
/> или />.
Чем меньше норма матрицы G, тем быстрее сходится итерационныйпроцесс.
Преобразование системыможно осуществить, просто решая каждое i-е уравнение относительноxi :
/>.
Метод Якоби использует следующий алгоритмпостроения приближений:
/>.
Если A — матрица с доминирующей диагональю, т.е. />, то метод Якоби сходитсяпри любом начальном приближении x(0).
Метод Якоби относится к одношаговымитерационным методам, когда для нахождения x(k+1) требуется помнить только однупредыдущую итерацию x(k). Для исследования сходимости удобнеезаписывать итерационные методы не в координатной, а в матричной форме,придерживаясь стандартной формы записи итерационных методов.
Канонической формойодношагового итерационного метода решения СЛАУ называется его запись в виде
/>,
где Bk+1 — матрица, задающая тот или иной итерационный метод, tk+1 — итерационный параметр. Числовыепараметры tk вводят для ускорения сходимости.Способ выбора итерационных параметров определяется при исследовании сходимостиметода, когда выясняется при каких значениях параметров метод сходится и когдасходимость будет наиболее быстрой (соответствующие параметры называютсяоптимальными).
Итерационный методназывают явным, если Bk+1 — единичная матрица. Неявныеитерационные методы имеет смысл применять лишь в том случае, когда решениесистемы уравнений с матрицей Bk требует меньше машиннойпамяти или времени или алгоритмически проще, чем решение исходной системы.
Методом простойитерации называютявный метод с постоянм параметром
/>, или/>,
где r(k) = Ax(k)-b — вектор невязки. Методсходится для симметричных положительно определенных матриц при />.
Дляокончанияитерационного процесса используют три способа. При первом определяютвеличину стабилизации и прекращают вычисления, если она меньше e, т.е.
/>.
Недостатком этого способаявляется то, что при медленно сходящихся итерациях величина стабилизации можетбыть малой, хотя приближенное решение сильно отличается от точного.
При втором способе вычисляютнормы невязки до начала итераций и на каждой итерации. Итерации прекращают привыполнении неравенства
/>.
При третьем способепредварительно оценивается число итераций, необходимое для получения заданнойточностиe.Если для погрешности итерационного метода выполняются оценки
/>,
где q (0,1), тометод сходится со скоростью геометрической прогрессии со знаменателемq.Можно определить, потребовав, чтобы qnn, достаточное для того, чтобыначальная погрешность уменьшилась в заданное число раз:
/>.
Целая часть числа n0(e) является минимальным числом итераций, необходимымдля получения заданной точности e.
Величина ln(1/q)является скоростью сходимости итерационного метода.
2. Описаниеиспользуемого метода
Для решения методомЗейделя система линейных алгебраических уравнений Ax= bдолжна быть приведена к виду x= Gx+f, где G — некоторая матрица, f-преобразованный вектор свободных членов. Затем выбирается начальное приближение- произвольный вектор x(0) — и строится рекуррентная последовательность векторов x(1), x(2),...,x(k),... по формуле
/>.
Для сходимости этойпоследовательности при любом начальном приближении необходимо и достаточно,чтобы все собственные значения матрицы G были по абсолютной величине меньше единицы. Напрактике это трудно проверить, и обычно пользуются достаточными условиями сходимости — итерации сходятся, если какая-нибудь норма матрицы меньше единицы, т.е.
/> или />.
Чем меньше норма матрицы G, тем быстрее сходится итерационныйпроцесс.
Преобразование системыможно осуществить, просто решая каждое i-е уравнение относительноxi :
/>.
Метод Зейделя использует следующий алгоритмпостроения приближений:
/>
Если A — матрица с доминирующей диагональю, т.е. />, то метод Зейделя сходитсяпри любом начальном приближении x(0).
/>
Метод Зейделя сходитсяпримерно так же, как геометрическая прогрессия со знаменателем ||G||. Если норма матрицы G близка к 1, то скорость сходимостиочень медленная. Для ускорения сходимости используется метод релаксации.Суть его в том, что полученное по методу Зейделя очередное значениепересчитывается по формуле:
Здесь 0нижняя релаксация,если w>1– верхняя релаксация. Параметр w подбирают так, чтобы сходимость метода достигалась заминимальное число итераций.
Метод Зейделя является одношаговымитерационным методам, когда для нахождения x(k+1) требуется помнить только однупредыдущую итерацию x(k).
Погрешность итерациивычисляется по формуле:
/>
n — порядок матрицы A.
Если d меньше заданной точности e, то итерационный процесс прекращают.
Элементы главнойдиагонали называются главными. Заметим, что если в ходе расчётов по данномуалгоритму на главной диагонали окажется нулевой элемент, то произойдет сбойпрограммы. Для того, чтобы избежать этого, следует перестановку строк такимобразом, чтобы на главной диагонали находились максимальные элементы строк. Т.е., если в k-й строке максимальным является i-йэлемент, необходимо поменять местами k-ю и i-ю строки, и поменять местамисоответствующие элементы вектора b.Такой выбор главного элемента необходим для сходимости итерационного процесса.
Приведёмблок-схему реализации данного метода:
/>
/>
/>
/>
/>
/>
/>
/>
/> /> /> /> /> /> /> /> /> /> /> /> /> /> /> /> /> />
/> /> /> /> /> /> /> /> />
/> /> /> /> /> /> /> /> />
3. Анализ результатов.
Скоростьсходимости итерационного процесса зависит от заданной матрицы коэффициентов. Взависимости от вида исходных данных( матрицы коэффициентов и матрицы b) программа подбирает оптимальный параметр релаксации w(при котором решениедостигается за минимальное число итераций).Длядостижения наивысшей скорости сходимости итерационного процесса для уравнения,заданного на рис.3 программой был выбран параметр релаксации w=1,26. Таким образом, былаприменена верхняя релаксация. Заданная точность e=0,0001 была достигнута за40 итераций. Графикзависимости количества итераций от параметра релаксации приведен на рис 1./> />
Рис. 1Длядостижения наивысшей скорости сходимости итерационного процесса для уравнения,заданного на рис.4 программой был выбран параметр релаксации w=0,98. Таким образом, былаприменена нижняя релаксация. Заданная точность e=0,0001 была достигнута за17 итераций. График зависимости количества итераций от параметра релаксацииприведен на рис 2.
/>
Рис. 2Правильностьрешения СЛАУ была проверена с помощью программного пакета Mathcad2000 professional. Отметим, что программа даётправильное решение СЛАУ почти во всех случаях, когда каждый элемент главнойдиагонали является максимальным в своей строке.
Вывод
Программа, разработаннаяв данной курсовой работе, реализует метод Зейделя для решения СЛАУ 6-гопорядка. Она даёт гарантированно правильное решение системы линейных уравнений,если каждый элемент главной диагонали матрицы коэффициентов являетсяединственным максимальным в своей строке, ненулевым, либо справедливы условия:максимальный элемент строки является единственным максимальным в своём столбце,ненулевым, а ни один из остальных элементов столбца не является максимальным всвоей строке, все элементы каждой строки кроме максимального одинаковы.
При исходных данных:
/>/>
была достигнута точность0,0001 в решении:/> />
за 2 итерации при параметререлаксации w=0,97.
Программа строит графикзависимости количества итераций от параметра релаксации для данной СЛАУ,находит параметр релаксации w, при котором решение достигается за минимальное количество итераций и,разумеется, само решение. Программа проста в эксплуатации и нетребовательна кресурсам. Реализованная в современной среде разработки Delphi 5.0, она без труда может быть доработана илиисправлена.
Недостатки программы: 1)применима не для всех систем линейных уравнений; 2)оптимальный параметррелаксации wвычисляется методом подбора, и, поэтому, количество итераций, требуемое для егоотыскания достаточно велико(около 18000), однако, для современных ПК, это неявляется затруднением.
Список использованнойлитературы
1. ВолковЕ.А. Численные методы. ¾ М.: Наука, 1987. ¾ 254 с.
2. КалиткинН.Н. Численные методы. ¾ М.: Наука, 1978. ¾ 512 с.
3. МудровА.Е. Численные методы для ПЭВМ на языках БЕЙСИК, ФОРТРАН и ПАСКАЛЬ. ¾ Томск, МП«Раско», 1992. ¾270 с.
4. СамарскийА.А., Гулин А.В. Численные методы. ¾ М.: Наука, 1989. ¾432с.
5. КэнтуМ. Delphi 4 для профессионалов ¾ СПб: «Питер», 1999 ¾1200с.
6. Delphi 5.0 help.
Приложение(распечаткапрограммы, результатов)
Распечаткапрограммы:
unit kurs1;
interface
uses
Windows,Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
Grids,StdCtrls, ComCtrls, ToolWin, Menus, Unit1, TeEngine, Series,
ExtCtrls,TeeProcs, Chart;
type
TFormk1 =class(TForm)
StringGrid1:TStringGrid;
Label2:TLabel;
Label3:TLabel;
Label4:TLabel;
StringGrid2:TStringGrid;
Label5:TLabel;
Label6:TLabel;
StringGrid3:TStringGrid;
Label7:TLabel;
Label8:TLabel;
Button1:TButton;
MainMenu1:TMainMenu;
Chart1:TChart;
N2: TMenuItem;
N3: TMenuItem;
N4: TMenuItem;
N1: TMenuItem;
Label1: TLabel;
Series1:TFastLineSeries;
procedureFormCreate(Sender: TObject);
procedurematrix;
procedureButton1Click(Sender: TObject);
procedureN1Click(Sender: TObject);
procedureN3Click(Sender: TObject);
procedureN4Click(Sender: TObject);
proceduredecision;
private
{ Privatedeclarations }
public
{ Publicdeclarations }
end;
var
Formk1:TFormk1;
// formk1:Tmainmenu;
implementation
var
n,m,i,j,k,l,number_of_iteration,min:integer;
delta,E,sum,max,W,tmp:extended;
A:array[1..6,1..6] of extended;
B:array[1..6] of extended;
X:array[1..6] of extended;
Xp:array[1..6] of extended;
am:array[1..200] of integer;
W_all:array[1..200]of extended;
procedureTFormk1.matrix;
begin
randomize;
for i:=1 to ndo stringgrid1.cells[i-1,0]:='*X'+inttostr(i);
for i:=0 ton-1 do
for j:=1 to mdo
StringGrid1.cells[i,j]:='2';
for i:=0 ton-1 do
StringGrid1.cells[i,i+1]:='3';
end;
{$R *.DFM}
procedureTformk1.decision;
begin
delta:=E+1;
number_of_iteration:=0;
for i:=1 to 6do X[i]:=B[i]/A[i,i];
while(delta>E) and (number_of_iteration
begin
for i:=1 to 6do Xp[i]:=X[i];
for i:=1 to 6do
begin
sum:=0;
for j:=1 to 6do sum:=sum+A[j,i]*X[j];
X[i]:=W*(B[i]-sum + A[i,i]*X[i])/A[i,i] + (1-W)*Xp[i];
end;
delta:=abs(X[1]-Xp[1]);
for i:=1 to 6do
if abs(X[i]-Xp[i])>deltathen delta:=abs(X[i]-Xp[i]);
inc(number_of_iteration);
end;
end;
procedureTFormk1.FormCreate(Sender: TObject);
begin
n:=6;m:=6;
matrix;
randomize;
stringgrid2.cells[0,0]:='*1';
for j:=1 to mdo
StringGrid2.cells[0,j]:='5';
end;
procedureTFormk1.Button1Click(Sender: TObject);
begin
series1.clear;
for i:=0 ton-1 do
for j:=1 to m do
A[i+1,j]:=strtofloat(StringGrid1.cells[i,j]);
for j:=1 to mdo
B[j] :=strtofloat(StringGrid2.cells[0,j]);
for i:=1 to 6do
begin
max:=abs(A[1,i]);
for j:=1 to 6do
ifabs(A[j,i])>=abs(max) then
begin
max:=A[j,i];
m:=j;
end;
if mithen
begin
for l:=1 to 6do
begin
tmp:=A[l,m];
A[l,m]:=A[l,i];
A[l,i]:=tmp;
end;
tmp:=b[m];
b[m]:=b[i];
b[i]:=tmp;
end;
end;
E:=0.0001;
W:=0.2;
l:=0;
whileW
begin
decision;
inc(l);
am[l]:=number_of_iteration;
W_all[l]:=W;
series1.addxy(W,number_of_iteration,'',clteecolor);
W:=W+0.01;
end;
min:=am[1];
for i:=1 to200 do
if(am[i]0) then
begin
min:=am[i];
W:=W_all[i];
end;
decision;
if(number_of_iteration>100) or (delta>E) then
begin
label2.Caption:='Программа не может решить данную СЛАУ.';
label3.Visible:=false;
end
else
begin
Chart1.BottomAxis.Automatic:=false;
Chart1.BottomAxis.minimum:=0.2;
Chart1.BottomAxis.maximum:=1.8;
Chart1.BottomAxis.increment:=0.1;
Chart1.LeftAxis.Automatic:=false;
Chart1.LeftAxis.minimum:=0;
Chart1.LeftAxis.maximum:=100;
Chart1.LeftAxis.increment:=5;
label6.visible:=false;
label7.visible:=true;
label8.visible:=true;
label1.visible:=true;
StringGrid3.visible:=true;
stringgrid3.cells[0,0]:='*1';
for i:=1 to 6do
StringGrid3.cells[0,i]:=floattostr(X[i]);
end;
end;
procedureTFormk1.N1Click(Sender: TObject);
begin
close;
end;
procedureTFormk1.N3Click(Sender: TObject);
begin
chart1.visible:=true;
end;
procedureTFormk1.N4Click(Sender: TObject);
begin
chart1.Visible:=false;
end;
end.
Результаты, рис. 3 и 4:
/>
Рис. 3
/> />
Рис. 4