Geum.ru

И. И. Мечникова лаборатория кафедра компьютерных методов экспериментальной экспериментальной физики физики компьютерный практикум

Вид материалаПрактикум

Содержание


Работа №7 Решение систем линейных алгебраических уравнений 1. Краткая теория
2. Реализация метода исключения Гаусса в библиотечном модуле Numerics
3. Примеры использования процедуры LinSys
Нулевой определитель системы
В а р и а н т 1
В а р и а н т 2
В а р и а н т 3
В а р и а н т 4
Подобный материал:


ОДЕССКИЙ ГОСУНИВЕРСИТЕТ им. И.И.Мечникова



ЛАБОРАТОРИЯ
КАФЕДРА


КОМПЬЮТЕРНЫХ МЕТОДОВ
ЭКСПЕРИМЕНТАЛЬНОЙ

ЭКСПЕРИМЕНТАЛЬНОЙ ФИЗИКИ
ФИЗИКИ



КОМПЬЮТЕРНЫЙ ПРАКТИКУМ


ДЛЯ ФИЗИКОВ


МЕТОДИЧЕСКИЕ УКАЗАНИЯ ДЛЯ СТУДЕНТОВ

ФИЗИЧЕСКОГО ФАКУЛЬТЕТА


Автор доцент П.А.Виктор


Р а б о т а 7


РЕШЕНИЕ СИСТЕМ ЛИНЕЙНЫХ

АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ


ОДЕССА - 1995




Работа №7

Решение систем линейных

алгебраических уравнений




1. Краткая теория



Системой n линейных алгебраических уравнений с n неизвестными называется система вида


µ § (1)


где µ § - неизвестные. Если определитель, составленный из коэффициентов системы,


µ §


не равен нулю, то система (1) имеет единственное решение. Это решение, в частности, может быть найдено с помощью правила Крамера:

µ §,

где µ § - определитель, получающийся из D при замене элементов k - го столбца µ § соответствующими свободными членами µ §. Однако формула Крамера удобна для теоретического исследования свойств решения, но весьма невыгодна для его численного нахождения. Дело в том, что для непосредственного раскрытия определителей системы с n неизвестными требуется порядка n!n арифметических операций, и уже при n=30 такое число операций недоступно современным компьютерам.

Одним из наиболее известных методов решения систем линейных алгебраических уравнений является метод исключения Гаусса. Прежде, чем рассматривать этот метод, вспомним два важных свойства систем линейных уравнений:

1. Решение системы не изменится, если правую и левую часть каждого уравнения умножить на любое число, отличное от нуля.

2. Решение системы не изменится, если к любому уравнению системы почленно прибавить (или вычесть) другое уравнение этой системы.

Метод исключения Гаусса основан на приведении системы к так называемому треугольному виду на основе использования упомянутых свойств. Рассмотрим его на примере следующей системы:


µ § (2)


Вычтем из второго уравнения системы первое, умноженное на такое число, чтобы уничтожился коэффициент при µ § (это число 1/2). Точно также добьемся исчезновения первого слагаемого в третьем уравнении, вычитая из него первое уравнение, умноженное на 7/4. В результате получим:

µ §


Аналогично в третьем уравнении избавимся от слагаемого с µ §, вычитая из него второе уравнение, умноженное на -17/10.


µ § (3)


Мы получили систему, приведенную к треугольному виду, решение которой совпадает с решением исходной системы (2). На этом завершается первая часть метода Гаусса, которая называется прямым ходом.

Вторая часть метода Гаусса - обратный ход - состоит в последовательном нахождении неизвестных, начиная с последнего уравнения системы (3). Так, в нашем случае из третьего уравнения непосредственно получаем µ §, подставляя это значение во второе уравнение, находим µ §, и, наконец, зная эти два неизвестных, из первого уравнения имеем µ §.

2. Реализация метода исключения Гаусса в библиотечном модуле Numerics



Метод исключения Гаусса реализован в процедуре LinSys, входящей в библиотечный модуль Numerics. Ниже приводится текст той части модуля, которая имеет отношение к указанной процедуре. Пользователи могут пропустить данный текст и ограничиться лишь рассмотрением порядка обращения к процедуре LinSys, который описан в следующем разделе.


unit Numerics;

interface


{$N+,E+}

const

ArraySize = (2*MaxInt) div SizeOf (extended);

type

ArrayType = array [1..ArraySize] of extended;


procedure LinSys { pешение системы N линейных

алгебpаических уpавнений

с N неизвестными методом Гаусса

с выбоpом ведущего элемента }


(N: integer; { число уpавнений }

var Matr; { pасшиpенная матpица системы }

var Sol; { pешение системы }

var Det: extended);{ опpеделитель системы }


{ ============================================ }


implementation


const

Eps: extended = 0; { нач. значение машинного эпсилон }


{ -------------------------------------------- }


procedure LinSys

(N: integer; var Matr; var Sol; var Det: extended);


var

A: ArrayType absolute Matr;

X: ArrayType absolute Sol;

k, l, m: integer;

C: extended;


function Ind (i, j: integer): integer; { внутpенняя }

{ вычисляет положение элемента в массиве-стpоке }

{ по значениям индексов i,j матpицы A }

begin

Ind:= Succ(N)*Pred(i) + Succ(j)

end; { Ind }


begin

{ вычисление машинного эпсилон пpи пеpвом вызове }

if Eps = 0 then

begin

Eps:= 1;

while 1+Eps/2 <> 1 do

Eps:= Eps/2;

Eps:= 20*Eps

end;


{ инициализация pешения и опpеделителя }

for k:= 1 to N do

X [k]:= 0;

Det:= 1;


{ пpямой ход }

for k:= 1 to N-1 do

begin

C:= Abs (A [Ind (k,k)]);

l:= k;

for m:= k+1 to N do

if Abs (A [Ind (m,k)]) > C then

begin

C:= Abs (A [Ind (m,k)]);

l:= m

end;


if C < Eps then

begin

Writeln ('НУЛЕВОЙ ОПРЕДЕЛИТЕЛЬ СИСТЕМЫ,');

Writeln ('ЕДИНСТВЕННОГО РЕШЕНИЯ НЕТ');

Halt

end;


if l <> k then { пеpестановка стpок }

begin

for m:= 0 to N do

begin

C:= A [Ind (k,m)];

A [Ind (k,m)]:= A [Ind (l,m)];

A [Ind (l,m)]:= C

end;

Det:= -Det

end;


for m:= k+1 to N do

begin

C:= A [Ind (m,k)] / A [Ind (k,k)];

A [Ind (m,k)]:= 0;

for l:= k+1 to N do

A [Ind (m,l)]:= A[Ind (m,l)]-C*A[Ind (k,l)];

A [Ind (m,0)]:= A[Ind (m,0)]-C*A[Ind (k,0)]

end

end;


{ обpатный ход }

for k:= N downto 1 do

begin

C:= 0;

for l:= k+1 to N do

C:= C + A [Ind (k,l)] * X [l];

X [k]:= (A [Ind (k,0)] - C) / A [Ind (k,k)]

end;


{ вычисление опpеделителя }

for k:= 1 to N do

Det:= Det * A [Ind (k,k)]


end; { LinSys }


END.

3. Примеры использования процедуры LinSys



Перед решением системы необходимо записать ее расширенную матрицу. Левый столбец этой матрицы (будем считать его нулевым) составляется из свободных членов µ §, а столбцы с 1-го по n-й образованы коэффициентами при неизвестных:1


µ §


Например, для системы (2) расширенная матрица записывается так:


µ §


Обращение к процедуре LinSys имеет вид:


LinSys (N, Matr, A, Det);

Здесь

N - число уравнений системы, integer;

Matr - расширенная матрица коэффициентов системы. Элементы матрицы должны иметь тип extended.

A - массив из N элементов типа extended. Сюда процедура LinSys записывает вычисленное решение системы;

Det - переменная типа extended. Сюда процедура записывает вычисленный определитель системы. В случае, если определитель системы равен нулю, процедура выводит на экран сообщение:

НУЛЕВОЙ ОПРЕДЕЛИТЕЛЬ СИСТЕМЫ

ЕДИНСТВЕННОГО РЕШЕНИЯ НЕТ

и выполнение программы прекращается.


Ниже приводится программа решения системы (2) с использованием процедуры LinSys, а также результаты ее работы. Программа начинается с директив компилятору, обеспечивающих подключение математического сопроцессора для обработки чисел в формате extended или его программную эмуляцию. Далее указывается, что программа использует библиотечный модуль Numerics. В разделе описания констант задается число уравнений NEqu и константа-массив A, содержащая расширенную матрицу системы размером 3 строки на 4 столбца. В разделе описания переменных описывается массив X из трех элементов для записи решения системы и переменная D для записи вычисленного определителя. Целочисленная переменная i используется в дальнейшем в качестве счетчика при выводе значений корней.

Основная часть программы предельно проста. Поскольку все исходные данные уже заданы в разделе описания констант, она сразу начинается с обращения к процедуре решения системы, после чего следует не требующий пояснений цикл вывода значений неизвестных и вывод значения определителя системы.

Отметим, что чаще всего для решения задачи не требуется вычисления определителя, однако в любом случае для работы процедуры LinSys под него нужно выделять переменную в разделе описаний.


program Lin1;


{$N+,E+}

uses Numerics;


const

NEqu = 3; { число уpавнений системы }


A: array [1..NEqu, 0..NEqu] of extended =

{ pасшиpенная матpица }

( ( 13, 4, 3, 1 ),

( -3, 2, -1, -1 ),

( 0, 7, 1, -3 ) ) ;


var

X: array [1..NEqu] of extended;

D: extended;

i: integer;


BEGIN


LinSys (NEqu, A, X, D);


{ вывод коpней и опpеделителя }

for i:= 1 to NEqu do

Writeln ('X[', i, '] =', X[i]:7:4);

Writeln ('Опpеделитель pавен', D:8:4)


END.


X[1] = 1.0000

X[2] = 2.0000

X[3] = 3.0000

Опpеделитель pавен 22.0000


Следующий пример демонстрирует использование процедуры LinSys для расчета разветвленных электрических цепей. Пусть имеется электрическая цепь из источников ЭДС и резисторов, изображенная на рисунке. Требуется рассчитать значения токов через каждый резистор.


µ §


µ §


Расчет сводится к решению системы линейных уравнений, которые составляются на основании правил Кирхгофа и содержат в качестве неизвестных значения токов.

Согласно правилам Кирхгофа:

1. Для любого узла цепи сумма токов втекающих в узел, равна сумме токов, вытекающих из узла.

2. Алгебраическая сумма падений напряжения на всех элементах любого замкнутого контура равна алгебраической сумме всех ЭДС, действующих в этом контуре.

Вначале стрелками отметим направления токов, которые мы принимаем за положительные. Пусть µ § - ток, протекающий через резистор µ §, и т.д. Первое правило Кирхгофа применим для узла A, а второе - для контуров ABCD и ADEF. В результате получаем три следующих уравнения:

µ §


Расширенная матрица системы имеет вид:


µ §


Ниже приводится программа, решающая систему линейных уравнений с данной матрицей, и результаты ее работы. Программа не нуждается в особых комментариях. Обратим внимание лишь на то, что начиная с версии 6.0 Turbo-Pascal допускает использование арифметических выражений при описании констант (это использовано при описании матрицы A). Отметим также, что поскольку значения ЭДС заданы в вольтах, а сопротивлений - в омах, решение системы представляет собой токи, выраженные в амперах.


program Lin2;


{$N+,E+}


uses Numerics;


const

NEqu = 3; { число уpавнений системы }


{ Значения сопpотивлений и ЭДС }

R1 = 30; R2 = 20; R3 = 40;

E1 = 80; E2 = 45;


A: array [1..NEqu, 0..NEqu] of extended =

{ pасшиpенная матpица }

( ( 0, 1, 1, -1 ),

( E1+E2, 0, R2, R3 ),

( E1, -R1, R2, 0 ) ) ;

var

I: array [1..NEqu] of extended;

D: extended;

k: integer;


BEGIN


LinSys (NEqu, A, I, D);


{ вывод токов чеpез pезистоpы }

for k:= 1 to NEqu do

Writeln ('I', k,' =', I[k]:7:4, ' ампеp');

Readln

END.


I1 =-0.8846 ампеp

I2 = 2.6731 ампеp

I3 = 1.7885 ампеp


Очевидно, отрицательное значение тока через первый резистор означает, что ток через него протекает в направлении, противоположном указанному стрелкой на схеме.


________________________




4. Задание



В соответствии с предложенным вариантом, рассчитайте значения токов через все резисторы изображенной на схеме электрической цепи, пользуясь процедурой LinSys. Результаты расчета выведите на экран.

В а р и а н т 1



µ §


µ §

В а р и а н т 2



µ §


µ §


В а р и а н т 3



µ §


µ §


В а р и а н т 4



µ §


µ §


В а р и а н т 5


µ §


µ §


В а р и а н т 6



µ §


µ §


* * *



1В математической литературе столбец свободных членов обычно записывается справа.



'