4Основы численных методов

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

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

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

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

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

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

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

87

новных численных методах и алгоритмах, уметь оценивать точность получаемого решения и потребные ресурсы ЭВМ.

4.1 Точность представления данных при решении прикладных задач

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

По этим причинам при работе с данными важной является оценка абсолютной погрешности их представления.

Пусть а — точное значение числа, а* — приближенное его значение, погрешностью ∆ числа а называется разность между приближенным и точным значением числа

∆ = а* - а

абсолютной погрешностью ∆а числа а называется абсолютная величина погрешности

∆a = |∆| == |а*- а|

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

δa = a *a *−a = ∆a *a

Обычно при описании погрешности представления данных используется форма

а= а* ± ∆а.

Погрешности можно разделить на три класса.

1.Неустранимая погрешность, возникающая за счет неточностей в математическом описании задачи и неточности задания исходных данных.

2.Погрешность метода. Она возникает за счет использования приближенного метода решения задачи.

3.Вычислительная погрешность. Она происходит при вычислениях за счет округления чисел.

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

Считается что, п первых значащих цифр числа верные, если абсолютная погрешность этого числа не превышает половины единицы разряда выраженного n-й значащей цифрой, считая слева направо.

Правило округления:

если первая из отброшенных цифр меньше 5, то оставшиеся цифры оставляют без изменения;

если первая из отброшенных цифр больше 5, то последняя цифра увеличивается на единицу;

если первая из отброшенных цифр равна 5 и после нее есть отличные от нуля цифры, то последняя цифра увеличивается на единицу;

если первая из отброшенных цифр равна 5 и после нее все остальные равны нулю, то последняя из оставшихся цифр остается без изменения.

Если положительное приближенное число имеет п верных знаков, то его отно-

сительная погрешность не превосходит 10 n - 1.

88

Необходимо иметь в виду, что при выполнении приближенных вычислений число значащих цифр в промежуточных результатах обычно не превышает числа верных цифр более чем на 2 единицы.

Следует помнить:

предельная абсолютная погрешность суммы или разности чисел равна сумме предельных абсолютных погрешностей слагаемых.

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

предельная относительная погрешность произведения (частного) двух отличных от нуля чисел не превосходит суммы относительных предельных погрешностей сомножителей (делимого и делителя).

4.2Некоторые часто используемые численные методы

Впрактике решения прикладных задач используется большое количество различных численных методов. Сущность большинства методов и их программная реализация подробно освещается в [1 – 5].

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

4.2.1Интерполяция.

Многие зависимости, используемые при решении прикладных задач, представляются в табулированном виде (таблицами). В таких таблицах непрерывная зависимость функции y=f(x) от аргумента x представляется рядом пар значений x и y, отражающих ряд точек на линии y=f(x). При использовании таблиц часто возникает необходимость в определении приближенного значения y на интервалах, расположенных между оговоренными точками. Другими словами решается задача нахождения значений функции для тех значений аргумента, которые находятся между значениями, имеющимися в таблице. Такой процесс называют интерполяцией.

Например, интерполируемая функция y=f(x) представлена в таблице набором пар значений x и y заданных в узлах интерполяции, и необходимо приближенно определить значение функции y* для заданного значения аргумента x*.

Аргумент - x	x1	x2	x3	…	xi-1	xi	xi+1	…	xn-1	xn
Функция - y	y1	y2	y3	…	yi-1	yi	yi+1	…	yn-1	yn

Линейная интерполяция на отрезке. Суть линейной интерполяции заключается в замене данной функции y=f(x) на отрезке [xi, xi+1] линейной функцией. При этом xi<x*<xi+1.

y* = y	i	+	yi+1 −yi	(x* −x )	.	(3.1)
				i
			xi+1 −xi

Процесс интерполяции осуществляется в следующей последовательности: сначала определяется интервал [xi, xi+1] внутри которого располагается заданное значение аргумента x*, а затем производится расчет по формуле (3.1). Величина y* определяется с некоторой погрешностью, которая тем больше чем выше значение второй производ-

89

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

Квадратичная интерполяция на отрезке. При квадратичной интерполяции предполагается, что заданная функциональная зависимость у=f (х) заменяется квадрат-

ной функцией у* = а(х*)2 + bx* + с на отрезке [xi, xi+2] при условии xi<x*<xi+2. Для нахождения коэффициентов а, b, с решают систему линейных уравнений:

yi = a xi2 +b xi + c;

yi+1 = a xi2+1 +b xi+1 + c; yi+2 = a xi2+2 +b xi+2 + c.

Откуда получаем:

	a =	( yi+2		− yi+1 ) ( xi+1 − xi ) −( yi+1 − yi ) ( xi+2						− xi+1 )	;
		( xi2+2		− xi2+1 ) ( xi+1 − xi ) −( xi2+1 − xi2 ) ( xi+2						− xi−1 )
		( y	i+1	− y	i	) − a ( x2	− x 2	)
	b =					i+1	i		;
				( xi+1 − xi )

c = yi − a xi2 −b xi

(3.2)

(3.3)

Определив а, b, с, тем самым определяют функцию у = ax2 + bх + + с, а по ней вычисляют промежуточное значение

y*=а(х*)2+bх*+с.

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

Полиномиальная интерполяция в области определения функции. Перспек-

тивным представляется использование единой интерполирующей функции во всей области определения табулированной функции. Чаще всего используется полиномиальная интерполяция. В этом случае для описания зависимости у = f (x) используется алгебраический многочлен (полином). Известно, что для любого набора точек (xi, yi) i=1,2,3…n существует единственный интерполяционный многочлен степени n-1

P	= a	0	+ a	1	x + a	2	x 2 +L+ a	n−1	x n−1	,	(3.4)
n−1
который проходит через все узлы интерполяции, т. е.
		Pn−1 ( xi ) = yi					(i =1,2,3,Kn).

Для нахождения коэффициентов а0, а1, а2,…аn-1 необходимо решить систему состоящую из n линейных алгебраических уравнений вида:

a	0	+ a x + a	2	x 2	+L+ a	n−1	x n−1	= y	i	(i =1,2,3,...n) .	(3.5)
		1

Следует отметить, что значение порядка полинома на единицу меньше количества заданных точек.

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

90

4.2.2Аппроксимация.

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

Аппроксимация табулированной функции полиномом. В рассматриваемом

случае функциональную зависимость y=f(x), заданную таблично (табулированную), приближенно отражают (аппроксимируют) полиномом, проходящим возможно ближе к точкам с координатами (xi, yi), но не требуют совпадения значений искомого полинома и табулированной функции в точках (xi, yi). При подобной аппроксимации чаще всего используется метод наименьших квадратов.

Рассмотрим в качестве аппроксимирующей функции полином степени k

f ( x,a) =a	+a x+a		x2 +La xk .	(3.6)
0	1	2	k

Будем минимизировать сумму квадратов рассогласований s значений заданной и аппроксимирующей функций во всех n точках (xi, yi)

s = ∑n [f ( x , a ) − yi ]2 .	(3.7)
i =1

Согласно теории необходимым условием минимума функции s является равенство нулю ее частных производных:

						∂f(x,a)		= 0	(j = 0 ,1 ,2 ,...,k ) .	(3.8)
						∂a j
Развернув (3.8) получаем систему уравнений для определения a0, a1, a2,…,ak
	∂s			= 2 ∑n	(a 0 + a 1 x i + a 2 x i2 +L + a k x ik − y i ) = 0;
	∂a 0
				i=1
	∂s			= 2 ∑n	(a 0 + a 1 x i		+ a 2 x i2		+L + a k x ik − y i x i ) x i = 0;	(3.9)
	∂a 1
				i=1
LLLLLLLLLL							LLLLLLLLLL LL ;
	∂s			= 2 ∑n	(a 0 + a 1 x i + a 2 x i2 +L + a k x ik − y i ) x ik = 0.
	∂a k			i=1

которая приводится к системе линейных алгебраических уравнений (СЛАУ) с матрицей Грамма вида:

n			n		n	n
a0 n + a1 ∑xi		+a2 ∑xi2		+L+ ak ∑xik =		∑yi;
i=1			i=1		i=1	i=1
n	n		n		n	n
a0 ∑xi +a1	∑xi2		+a2 ∑xi3		+L+ak ∑xik+1	= ∑yixi;	(3.10)
i=1	i=1		i=1		i=1	i=1

LLLLLLLLLLLLLLLLLLLLLLLL;

n	n	n
a0 ∑xik +a1	∑xik+1 +a2	∑xik+2
i=1	i=1	i=1

n	n
+L+ak ∑xi2k	= ∑yixik.
i=1	i=1

Таким образом, для определения коэффициентов как интерполирующего, так и аппроксимирующего полиномов необходимо решить систему линейных алгебраических уравнений. Следует отметить, что многие математические модели тем или иным способом приводятся к СЛАУ и для пользователя ЭВМ полезно ознакомиться с основными методами решения СЛАУ.

91

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

Пусть А - квадратная матрица порядка п и b – вектор правых частей, состоящий из п компонент. Задача состоит в решении системы п уравнений:

Ах=b.

(3.11)

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

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

Запишем систему линейных алгебраических уравнений (3.11) в развернутом виде

a11 x1 + a12 x12 + a13 x3 +L+ a1n x n + = b1 ;
a 21 x1 + a 22 x12 + a 23 x3 +L+ a 2 n x n + = b2	;
LLLLLLLLLLLLLLLLLL;	(3.12)

a n1 x1 + a n 2 x12 + a n 3 x3 +L+ a nn x n + = bn

Необходимым и достаточным условием существования единственного решения такой системы является неравенство нулю ее определителя.

Наиболее распространенным из прямых методов решения таких уравнений является метод Гаусса и его модификации.

Метод исключения Гаусса. Как и многие другие, этот метод основан на приведении исходной системы к "треугольному" виду путем последовательного исключения неизвестных во входящих в эту систему уравнениях. Процесс исключения начинается с того, что все члены первого уравнения делят на а11 (эта процедура называется нормированием уравнения). Затем первое уравнение последовательно умножается на первые коэффициенты остальных уравнений системы, т. е. на коэффициенты аi1, и вычитается из каждого из этих (соответствующих) уравнений. По окончании этого процесса первая переменная (х1) будет исключена из всех уравнений, кроме первого уравнения.

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

менная (xn).

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

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

Количество арифметических операций в методе Гаусса связано с размерностью системы и примерно равно 2/3 п3.

92

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

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

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

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

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

Итерационные методы применяются к системам уравнений (3/12), предварительно преобразованным в форму

	x i =	1	( b i − ∑n	a ij x j )	i = 1,2,3,...n .	(3.13)
		a ii
			j =1
			j ≠ i

Итерации продолжаются до тех пор, пока для каждого неизвестного не выполнится

условие
		xm −xm−1			≤ε	,	(3.14)
где ε- абсолютный итерационный		i	i
		допуск,			m - номер итерации.

По способу организации итерационного уточнения значений искомых величин различают метод простых итераций (метод Якоби) и метод Зейделя (Гаусса-Зейделя).

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

	зультатам предыдущей итерации xm−1				,xm−1	,..xm−1 .		В этом случае счетная формула при-
	обретает вид:			1	2		n
	x im =	1	( b i − ∑n	a ij	x mj		− 1 )	i = 1,2,3,...n	(3.15)
		a ii
			j =	1
			j ≠	i
							93