Метод Зейделя. В методе Зейделя для вычисления на данной итерации значений каждой из искомых величин xim используются как значения, взятые по результатам

предыдущей итерации, так и уже уточненные в процессе данной итерации. Расчет ведется по формуле

		1	i −1	n
	xim =		( bi − ∑ a ij x mj	− ∑ a ij x mj	−1 ) i = 1,2,3,...n	(3.16)
		a ii
			j =1	j =i +1

Для абсолютной сходимости итерационных методов необходимо, чтобы значение каждого из диагональных элементов матрицы СЛАУ были преобладающими по абсолютной величине по сравнению с другими элементами строки. Условие сходимости можно обеспечить преобразованием исходной матрицы путем перестановки уравнений и неизвестных. При использовании любого итерационного метода итерационный процесс завершается, когда выполняются условия (3.14) для всех искомых значений.

Методы простых итераций и Зейделя имеют разные области сходимости. В ряде случаев метод Зейделя обеспечивает более быструю сходимость итерационного процесса, нежели метод простых итераций.

Следует отметить, что величина абсолютного итерационного допуска ε не является погрешностью решения, но ее величина близка к достигаемой точности решения.

Контроль полученных решений можно провести путем подстановки их результатов в исходную СЛАУ и вычисления невязок s (разностей между правыми и левыми частями уравнений).

Невязка каждой i-ой строки определятся как

n
si =bi −∑aij xj	(3.17)

j=1

4.2.4 Решение систем нелинейных уравнений.

Рассмотренные методы можно применять и к итерационному решению систем нелинейных уравнений. При этом самым простым приемом является сведение системы нелинейных уравнений к системе квазилинейных уравнений, т.е. таких уравнений, которые по структуре соответствуют СЛАУ, но коэффициенты которых не постоянны, а определяются на каждой текущей итерации по результатам предыдущей. Другими словами в этом случае учет всех нелинейных зависимостей осуществляется путем периодического пересчета значений коэффициентов матрицы. Более подробно этот прием рассмотрен далее в примере № 8.

4.2.5 Решение алгебраических и трансцендентных уравнений.

Во многих научных и инженерных задачах часто возникает необходимость решения уравнений вида

f ( x ) = 0

(3.18)

где f - заданная функция.

Уравнение (3.18) может иметь как конечное, так и бесконечное количество решений, что зависит от вида функции f(x).

Решениями или корнями уравнения (3.18) называются такие значения х, которые при подстановке в уравнение (3.18) обращают его в тождество. Нахождение корней таких уравнений аналитическими методами (т.е. путем записи формулы, выражающей искомую величину х в явном виде) возможно только для простейших случаев. Поэтому чаще всего приходится решать уравнения вида (3.18) численными методами. Кроме того, иногда, даже при наличии аналитического решения, имеющего сложный вид, быва-

94

ет проще провести численное решение по известному алгоритму, чем программировать громоздкую аналитическую формулу.

Численное решение уравнения (3.18) обычно проводят в два этапа. На первом этапе необходимо отделить корни уравнения, т.е. найти такие интервалы изменения переменной х, где расположен только один корень. На втором этапе проводят уточнение отделенных корней, т.е. находят корни с заданной точностью. Для этого разработан богатый набор алгоритмов и программ.

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

Для уточнения значений отделенных корней преимущественно используются методы половинного деления (дихотомии), хорд, простых итераций и Ньютона.

Метод половинного деления. Считаем, что отделение корней уравнения (3.18) проведено и на отрезке [а,b] расположен один корень, который необходимо уточнить с абсолютным итерационным допуском ε. Метод половинного деления, заключается в

следующем.

Как видно из рисунка 3.1,

вначале определяем середину от-

f(x)

резка [а, b] x*=(а+b)/2 и вычисляем

f(b)

функцию f(x*). Далее делаем выбор,

какую из двух частей отрезка взять

f(x*)

для дальнейшего уточнения корня.

Если левая часть уравнения f(x) есть

2ε

непрерывная функция аргумента х,

то корень будет находиться в той

2ε1

x

половине отрезка, на концах кото-

b

a

x*

рой f(x) имеет разные знаки. На ри-

корень уравнения

f(a)

сунке 3.1 это будет отрезок [а, х*].

Для очередного шага уточнения

Рис. 3.1. Метод половинного деления

точку b перемещаем в точку х* и

продолжаем процесс аналогичным

образом. Итерационный процесс будем продолжать до тех пор, пока интервал [а, b] не станет меньше заданного итерационного допуска ε.

Используя метод половинного деления, следует учитывать, что функция f(x) вычисляется с некоторой абсолютной погрешностью, а вблизи корня значения функции f(x) малы по абсолютной величине и поэтому могут оказаться сравнимыми с погрешностью ее вычисления. Другими словами, при подходе к корню существует опасность попасть в полосу “шумов” 2ε1 (рис. 3.1) и дальнейшее уточнение корня окажется невозможным. В этом случае процесс уточнения значения корня прекращается. При практических расчетах обычно принимают ε1=ε.

Метод половинного деления обладает быстрой сходимостью, так как за каждую итерацию интервал, где расположен корень, уменьшается в два раза. Таким образом через п итераций интервал будет равен (b - a)/2n.

Метод хорд. Метод хорд также, как и метод половинного деления, предназначен для уточнения корня на интервале [a,b], на концах которого значения функции f(x) разного знака. Первое приближение берется в точке x1, где пересекается прямая линия, проведенная через точки f(a) и f(b) с осью x. В качестве нового интервала для продолжения итерационного процесса выбирается тот из двух [a,x1] или [x1,b], на концах которого функция f(x) имеет разные знаки. В нашем случае это интервал [a,x1]. На

95

следующих итерациях после уточнения
границ нового интервала расчет прово-	f(x)		f(b)
дится по формуле:	Корень уравнения
x=(f(b)*a - f(a)*b)/(f(b) - f(a));		f(x1)
Итерационный процесс уточнения	2ε
корня заканчивается, когда выполнится
условие |xi-xi-1|<ε или когда значение
функции f(xi) попадет в область шума	2ε1
|f(xi)|<ε1.	a	x1		b	x
Метод Ньютона (метод каса-		x2
тельных). Рассмотрим графическую	f(a)
иллюстрацию метода на рис. 3.3.
Предположим, что известно на-	Рис. 3.2 Метод хорд
чальное приближение x0 к корню. В
точке x0 вычислим левую часть решаемого уравнения f(xQ), а также производную в этой
точке f’’(х0). Следующее приближение к корню найдем в точке х1, где касательная к
функции f(x), проведенная из точки (x0,f0), пресекает ось абсцисс. Для этого используем
формулу (3.19) построенную на основе элементарных геометрических соображений.

x1

= x0 −

f ( x0 )

(3.19)

f ' ( x0

)

Затем считаем точку х1 в ка-

f(x)

честве начальной и продолжаем

f(x0)

итерационный процесс (определяем

следующее приближение x2)

и т.д.

корень уравнения

Из рис. 3.3 видно, что таким спосо-

бом можно быстро приблизиться к

f(x1)

корню х*. При этом с каждой ите-

рацией расстояние между очеред-

f(x2)

ным xm и предыдущим хm-1 прибли-

жениями к корню будет умень-

шаться. Расчет уточненного значе-

x2

x1

x0

x

x*

ния корня в общем случае прово-

Рис. 3.3 Метод Ньютона

дим по уравнению (3.20). Процесс

уточнения корня закончим, когда выполнится условие (3.21).

xm = xm−1 −

f ( xm−1

)

(3.20)

f ' ( xm−1

)

xm+1 - хm < ε,

(3.21)

где ε – абсолютный итерационный допуск, m – номер итерации (приближения).

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

Метод простых итераций. В этом методе решается уравнение (3.22).

x=φ(x)

(3.22)

Процесс решения заключается в следующем. Пусть имеется приближенное значение корня xm-1. Для получения уточненного значения используем итерационную формулу:

96

xm = φ(xm-1)

(3.23)

Уточнения (итерации) выполняются до тех пор, пока не выполнится условие (3.19). Процесс уточнений может быть сходящимся, когда от итерации к итерации непрерывно уменьшается значение xm+1 - хm или может быть расходящийся когда значение xm+1-хm непрерывно растет. Условием сходимости итерационного процесса является

φ′(x) < 1

(3.24)

Переход от уравнения (3.18) к уравнению (3.22) можно осуществить различными способами в зависимости от вида функции f(x). При таком переходе необходимо построить функцию φ(x) так, чтобы выполнялось условие сходимости (3.24). Чаще всего используется преобразование (3.25)

φ(x)= x+p f(x).

(3.25)

Преобразование (3.25) удобно тем, что, изменяя значение параметра p можно управлять скоростью сходимости итерационного процесса.

4.3 Использование подпрограмм при решении прикладных задач

Как отмечалось ранее, для часто используемых вычислительных методов разрабатываются подпрограммы (процедуры и функции) которые хранятся в библиотеках под уникальными именами. Из программы пользователя к библиотечным подпрограммам можно обратиться по их именам. В системе программирования Borland Pascal такого рода библиотеки называются модулями. Результатом компиляции исходного текста модуля является файл с расширением tpu (tpp, tpw), который хранится на жестком диске. Расширение зависит от платформы разрабатываемого приложения. Расширение tpu соответствует обычному режиму real mode application. Пользователь при выполнении своей программы обращается к модулю как к библиотеке.

Для использования подпрограмм пользователю в начале своей программы (после заголовочного оператора) необходимо указать имена модулей, содержащих используемые в данной программе процедуры и функции. Для этого записывается ключевое слово uses и перечень имен модулей, разделенных запятыми с символом; в конце списка. В тексте программы обращение к подпрограммам оговоренных модулей производится по именам процедур или функций.

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

Список формальных параметров размещается в круглых скобках после имени подпрограммы, помещенной в модуль. В теле этой подпрограммы имена формальных параметров используются для описания алгоритма вычисления. Для обращения к выбранной подпрограмме из программы пользователя используется имя подпрограммы после которого размещается (в скобках) список фактических параметров. Списки фактических и формальных параметров должны быть согласованы по количеству, порядку следования и типам параметров. Другими словами в подпрограмме формальные параметры используются как «псевдонимы» фактических параметров для описания сути вычислений, а при обращении к этой процедуре формальные параметры получают конкретные значения соответствующим им фактических параметров.

Формальные параметры подразделяются на параметры-значения и параметрыпеременные.

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

97

программы никак не влияет на значения соответствующих фактических параметров. Поэтому соответствующие им фактические параметры могут быть как переменными, так и константами.

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

В списках формальных параметров параметры-значения и параметры-переменные описываются разными способами. После списка однотипных параметров-значений ставится символ : и ключевым словом указывается тип параметров-значений. Перед списком параметров-переменных, кроме того, обязательно ставится ключевое слово var.

Для иллюстрации приемов работы с подобными библиотеками рассмотрим модуль Library1, содержащий вычислительные процедуры и функции, отражающие рассмотренные ранее некоторые численные методы. Исходный текст этого модуля (содержимое файла Librery1.pas) приводится в приложении.

Исходный текст модуля включает два основных раздела:

interface,

implementation.

В разделе interface приводится описание типов отдельных объектов (оно предваряется ключевым словом type) и перечень заголовочных операторов процедур и функций, помещенных в данный модуль.

Указанные описания типов необходимо использовать при разработке программы пользователя!

В разделе implementation размещаются тексты процедур и функций.

Чаще всего для использования предопределенных процедур пользователю достаточно ознакомиться с разделом interface.

Рассмотрим раздел interface модуля Library1. unit Library1;

interface

type

Mat = array[1..20,1..21] of Real; Vec = array[1..101] of Real;

Xfn = function(X:Real):Real;

procedure Gauss(N:Integer; A:Mat; var X:Vec; var S:Real); procedure Zeidel(N:Integer; A:Mat; var X:Vec; var E:Real;

var M:Integer);

procedure VvodA(var N:Integer; var A:Mat); procedure Testing(N:Integer; A:Mat; X:Vec); procedure GRAM(N:Integer; var A:Mat; var X,Y:Vec;

var M,K:Integer);

procedure Newton(var X:Real; E:Real; F_div_df:Xfn); procedure Dix(A,B,E:Real; var X:Real; F:Xfn; var Kd:Integer); procedure Chord(A,B,E:Real; var X:Real; F:Xfn;

98

procedure	var Kd:Integer);
	krest(x,y:longint);
procedure	Grafic(nn:integer; var xx:vec; var yy:vec;
function	n:integer; var x:vec; var y:vec);
	Power(Xx:Real; P:Integer):Real;
function	Deg_Rad(Gr:Real):Real;

Процедура Gauss предназначена для решения СЛАУ методом Гаусса с выбором главного элемента. Формальные параметры: N – порядок матрицы, A – матрица размером (N,N+1), X – вектор решения, S – признак вырожденности матрицы (если S=0

–матрица вырождена).

Врассматриваемой матрице последний столбец (A[1,N+1] A[N,N+1]) включает вектор правых частей. Такой прием используется для упрощения обмена данными.

Процедура Zeidel предназначена для решения СЛАУ итерационным методом

Зейделя. Формальные параметры: N – порядок матрицы, A – матрица размером (N,N+1), X – вектор решения, E – абсолютный итерационный допуск, M – количество выполненных итераций.

Процедура VvodA предназначена для ввода значений коэффициентов СЛАУ и вектора правых частей. Формальные параметры: N – порядок матрицы, A – матрица размером (N,N+1).

Процедура Testing предназначена для проверки решения СЛАУ путем сравнения заданного вектора правых частей и полученного при подстановке решения в СЛАУ. Формальные параметры: N – порядок матрицы, A – матрица размером (N,N+1), X – вектор решения.

Процедура Gram предназначена для формирования матрицы Грамма. Формальные параметры: N – количество точек табулированной функции, A – матрица Грамма размером (M,M), X,Y – массивы координат точек табулированной функции, M – порядок матрицы Грамма, K – порядок аппроксимирующего полинома.

Процедура Newton предназначена для итерационного уточнения корня уравнения методом Ньютона. Формальные параметры: X – значение корня, E – абсолютный итерационный допуск, F_div_df – параметр соответствующий отношению функции к ее производной.

Процедура Dix предназначена для итерационного уточнения корня уравнения методом половинного деления. Формальные параметры: A,B – границы интервала на котором ищется уточненное значение корня уравнения, E – абсолютный итерационный допуск, X – значение корня, F – функция, Kd – признак удачного завершения итерационного процесса.

Процедура Chord предназначена для итерационного уточнения корня уравнения методом хорд. Формальные параметры: A,B – границы интервала на котором ищется уточненное значение корня уравнения, E – абсолютный итерационный допуск, X – значение корня, F – функция, Kd – признак удачного завершения итерационного процесса.

Процедура Krest предназначена для вывода на графике крестообразного маркера точки. Используется в процедуре Grafic . Формальные параметры: X,Y – координаты центра маркера.

Процедура Grafic предназначена для построение двух графиков однопараметрических зависимостей. Один график отражается крестообразными маркерами, другой ломаной линией. Формальные параметры: Nn – количество точек, отражаемых маркерами, Xx, Yy – массивы координат точек отражаемых маркерами, N - количество точек излома линии графика, X, Y – массивы координат точек излома линии графика. Для удобства пользователя на графике значения Y=0 и X=0 отражаются жирными линиями.

99