Роман Гордон
дата публикации 03-07-2008 11:52

---

Версия для печати

1. Преамбула

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

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

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

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

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

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

2. Сущность МНК (с примером).

Сущность обоснования МНК (по Гауссу) заключается в допущении, что «убыток» от замены точного (неизвестного) значения величины ее приближенным вычисленным значением пропорционален квадрату ошибки. Тогда оптимальной оценкой естественно считать такую величину, для которой среднее значение «убытка» минимально.

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

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

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

Известно, что любому измерению присуща случайная погрешность. Из графика снятой зависимости параметра от времени видно, что она приблизительно линейная. Требуется найти прямую, которая поможет выяснить значение искомого параметра в произвольный момент времени. Уравнение прямой y = a·x + b, где в нашем случае х – секунда отсчета (совпадающая с номером отсчета), y – значение параметра на этом отсчете. Чтобы найти оптимально подходящую прямую, нужно найти такие коэффициенты a и b, при которых сумма ∑([a·x_i + b − y_i]²) по всем отсчетам i была бы минимальной. Здесь a·x_i + b задает неизвестное значение прямой на отсчете i, y_i – известное измеренное значение, выражение под суммой – квадрат отклонения измеренной величины от искомой оценочной. Итак, имеем функцию

Q(a, b) = ∑([a·x_i + b − y_i]²),

и нужно найти такие значения a и b, при которых функция имеет минимум. Известно, что в точке минимума первая производная функции равна нулю. У нас две переменные, значит

dQ/da = ∑(2·[a·x_i + b − y_i]·x_i)
dQ/db = ∑(2·[a·x_i + b − y_i]·1).

Раскроем скобки, приравняем производные к нулю и получим систему уравнений с двумя неизвестными а и b:

a·∑(x_i²) + b·∑(x_i) = ∑(x_i·y_i)
a·∑(x_i) + b·∑(1) = ∑(y_i)

Понятно, что ∑(1) равна количеству отсчетов. Для наглядности попробуем задать серию из десятка точек «с потолка» с приблизительно линейной зависимостью y от x:

x	1	2	3	4	5	6	7	8	9	10
y	1	2	2	3	3	3	4	3	5	5

Тогда система выглядит так:

385·a + 55·b = 203
55·a + 10·b = 31

отсюда b = (31 − 55·a) / 10 = 3,1 − 5,5·a, тогда 385·a + 170,5 − 302,5·a = 203, значит, a = 0,(39), b = 0,9(3). Результат можно увидеть на рисунке 1(а). Для еще большей наглядности проверим более явную линейную зависимость, заданную таблично:

x	1	2	3	4	5	6	7	8	9	10
y	1	1	2	2	3	3	4	4	5	5

Здесь система получается такая:

385·a + 55·b = 205
55·a + 10·b = 30

отсюда а = 0,(48), b = 0,(3), и можно посмотреть, как это выглядит на рисунке 1(б).

a) б)

Рисунок 1. Пример аппроксимации табличных точек прямой

3. Описание аппроксимации плоскостью.

Имеется поверхность, заданная точками в виде двумерного массива. То есть, индексы массива представляют собой координаты x и y, а значения элементов массива – значения координаты z соответствующих точек.

Необходимо найти такую функцию a·x + b·y + c = z, которая будет оптимальной с точки зрения МНК, то есть вычислить соответствующие коэффициенты a, b и c. Разность между аппроксимированным и измеренным значениями в конкретной точке [x_i,j, y_i,j] равна [a·x_i,j + b·y_i,j + c − z_i,j], где x_i,j и y_i,j – координаты точки (в нашем случае x_i,j = i, а y_i,j = j), a·x_i,j + b·y_i,j + c – значение z на искомой плоскости в точке [x_i,j, y_i,j], а z_i,j – известное измеренное значение (элемент массива). Нужно найти такие a, b, c, при которых функция

Q(a, b, c) = ∑([a·x + b·y + c − z]²)

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

Подразумевается сумма по всем элементам двумерного массива, то есть

Q(a, b, c) = ∑_i∑_j([a·x_i,j + b·y_i,j + c − z_i,j]²),

где i изменяется от 0 до количества столбцов массива, j – от 0 до количества строк, х_i,j совпадает с i, y_i,j совпадает с j. Для простоты будет использоваться сокращенная запись. Итак, производные по a, b и с:

dQ/da = ∑(2·[a·x + b·y + c − z]·x)
dQ/db = ∑(2·[a·x + b·y + c − z]·y)
dQ/dc = ∑(2·[a·x + b·y + c − z]·1),

значит, система уравнений будет выглядеть так:

a·∑(x2)	+ b·∑(x·y)	+ c·∑(x)	= ∑(x·z)
a·∑(x·y)	+ b·∑(y2)	+ c·∑(y)	= ∑(y·z)
a·∑(x)	+ b·∑(y)	+ c·∑(1)	= ∑(z)

Очевидно, что ∑(1) равна числу элементов массива.

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

1. перемена местами двух уравнений в системе;
2. умножение какого-либо уравнения системы на действительное число, не равное 0;
3. прибавление к одному уравнению другого уравнения, умноженного на произвольное число.

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

Приступим к преобразованию системы уравнений к треугольному виду. А для удобства перепишем ee несколько иначе. a, b и с – неизвестные, и традиционно обозначим их через х, то есть х1, х2, х3. Множители, выраженные через суммы, фактически известны (могут быть получены из имеющихся данных путем элементарных вычислений). Обозначим их буквами a, b и c c индексами, указывающими на номер уравнения в системе. Правая часть не содержит неизвестных, обозначим ее буквой d, тоже с индексом:

a1 = ∑(x2)	b1 = ∑(x·y)	c1 = ∑(x)	d1 = ∑(x·z)	для первого уравнения
a2 = b1 = ∑(x·y)	b2 = ∑(y2)	c2 = ∑(y)	d2 = ∑(y·z)	для второго уравнения
a3 = с1 = ∑(x)	b3 = с2 = ∑(y)	c3 = ∑(1)	d3 = ∑(z)	для третьего уравнения

Таким образом, система примет вид:

a1·x1 + b1·x2 + c1·x3 = d1
a2·x1 + b2·x2 + c2·x3 = d2
a3·x1 + b3·x2 + c3·x3 = d3

Избавимся от х1 во втором и третьем уравнениях системы, воспользовавшись третьим эквивалентным преобразованием. Определим такие коэффициенты:

k12 = −a2/a1 – коэффициент для преобразования второго уравнения, k13 = −a3/a1 – коэффициент для преобразования третьего.

Теперь если умножить каждый элемент первого уравнения на k12 и получившеся уравнение прибавить ко второму, множитель при х1 второго уравнения обнуляется:

a1·k12·x1 + a2·x1 = (a1·(−a2/a1)·x1) + a2·x1 = (−a2·a1/a1 + a2)·x1 = 0·x1 = 0
b1·k12·x2 + b2·x2 = (b1·k12 + b2)·x2

Аналогично с k13 – умножим на него каждый член первого уравнения и прибавим первое уравнение к третьему. Система становится такой:

a1·x1 + b1·x2 + c1·x3 = d1
(b2 + k12·b1)·x2 + (c2 + k12·c1)·x3 = d2 + k12·d1
(b3 + k13·b1)·x2 + (c3 + k13·c1)·x3 = d3 + k13·d1

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

k23 = −(b3 + k13·b1) / (b2 + k12·b1)

Получаем систему, готовую к решению:

a1·x1 + b1·x2 + c1·x3 = d1
(b2 + k12·b1)·x2 + (c2 + k12·c1)·x3 = d2 + k12·d1
((c3 + k13·c1) + k23·(c2 + k12·c1))·x3 = (d3 + k13·d1) + k23·(d2 + k12·d1)

Решение очевидно:

x3 = ((d3 + k13·d1) + k23·(d2 + k12·d1)) / ((c3 + k13·c1) + k23·(c2 + k12·c1))
x2 = (d2 + k12·d1 − (c2 + k12·c1)·x3) / (b2 + k12·b1)
x1 = (d1 − b1·x2 − c1·x3) / a1

Подставив найденные х1, х2 и х3 вместо а, b и с в формулу a·x + b·y + c = z, получаем уравнение плоскости, аппроксимирующей нашу поверхность оптимально с точки зрения МНК.

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

4. Описание аппроксимации фигурой вращения

В качестве аппроксимирующей возьмем поверхность, задаваемую уравнением

z = e^{(a·x2 + b·y2 + c·x·y + d·x + f·y + g)}

Она имеет форму, напоминающую колокол. Для удобства прологарифмируем обе части уравнения, так удобнее будет дифференцировать:

ln(z) = a·x² + b·y² + c·x·y + d·x + f·y + g

Теперь получим функцию, представляющую собой сумму квадратов отклонений известных точек от соответствующих точек на искомой поверхности:

Q = ∑([a·x² + b·y² + c·x·y + d·x + f·y + g − ln(z)]²)

Опять же, необходимо найти такие a, b, c, d, f и g, при которых функция минимизирована, то есть ее частные производные

dQ/da = ∑(2·[a·x² + b·y² + c·x·y + d·x + f·y + g − ln(z)])·x²
dQ/db = ∑(2·[a·x² + b·y² + c·x·y + d·x + f·y + g − ln(z)])·y²
dQ/dc = ∑(2·[a·x² + b·y² + c·x·y + d·x + f·y + g − ln(z)])·x·y
dQ/dd = ∑(2·[a·x² + b·y² + c·x·y + d·x + f·y + g − ln(z)])·x
dQ/df = ∑(2·[a·x² + b·y² + c·x·y + d·x + f·y + g − ln(z)])·y
dQ/dg = ∑(2·[a·x² + b·y² + c·x·y + d·x + f·y + g − ln(z)])·1

равны нулю. После раскрытия скобок и приравнивания к нулю получим следующую систему уравнений:

a·∑(x4)	+ b·∑(x2·y2)	+ c·∑(x3·y)	+ d·∑(x3)	+ f·∑(x2·y)	+ g·∑(x2)	= ∑(x2·ln(z))
a·∑(x2·y2)	+ b·∑(y4)	+ c·∑(x·y3)	+ d·∑(x·y2)	+ f·∑(y3)	+ g·∑(y2)	= ∑(y2·ln(z))
a·∑(x3·y)	+ b·∑(x·y3)	+ c·∑(x2·y2)	+ d·∑(x2·y)	+ f·∑(x·y2)	+ g·∑(x·y)	= ∑(x·y·ln(z))
a·∑(x3)	+ b·∑(x·y2)	+ c·∑(x2·y)	+ d·∑(x2)	+ f·∑(x·y)	+ g·∑(x)	= ∑(x·ln(z))
a·∑(x2·y)	+ b·∑(y3)	+ c·∑(x·y2)	+ d·∑(x·y)	+ f·∑(y2)	+ g·∑(y)	= ∑(y·ln(z))
a·∑(x2)	+ b·∑(y2)	+ c·∑(x·y)	+ d·∑(x)	+ f·∑(y)	+ g·∑(1)	= ∑(ln(z))

Для удобства решения методом Гаусса, как и в предыдущем варианте с плоскостью, перепишем систему несколько иначе: a, b, c, d, f и g заменим на х1, х2, х3, х4, х5 и х6, а множители при них обозначим через а с индексом, обозначающим позицию в уравнении – строку и столбец. То есть, а11 = ∑(x⁴), а12 = ∑(x²·y²), ..., а21 = ∑(x²·y²), а22 = ∑(y⁴), ..., а66 = ∑(1). Правую часть системы заменим на b аналогичным образом : b1 = ∑(x²·ln(z)), ..., b6 = ∑(ln(z)).

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

a11·x1 + a12·x2 + a13·x3 + a14·x4 + a15·x5 + a16·x6 = b1
a21·x1 + a22·x2 + a23·x3 + a24·x4 + a25·x5 + a26·x6 = b2
a31·x1 + a32·x2 + a33·x3 + a34·x4 + a35·x5 + a36·x6 = b3
a41·x1 + a42·x2 + a43·x3 + a44·x4 + a45·x5 + a46·x6 = b4
a51·x1 + a52·x2 + a53·x3 + a54·x4 + a55·x5 + a56·x6 = b5
a61·x1 + a62·x2 + a63·x3 + a64·x4 + a65·x5 + a66·x6 = b6

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

k12 = −a21/a11

k13 = −a31/a11

k14 = −a41/a11

k15 = −a51/a11

k16 = −a61/a11

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

a2_22 = a22 + k12·a12	a2_23 = a23 + k12·a13	a2_24 = a24 + k12·a14	...	b2_2 = b2 + k12·b1
a2_32 = a32 + k13·a12	a2_33 = a33 + k13·a13	a2_34 = a34 + k13·a14	...	b2_3 = b3 + k13·b1
...

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

a11·x1 + a12·x2 + a13·x3 + a14·x4 + a15·x5 + a16·x6 = b1
a2_22·x2 + a2_23·x3 + a2_24·x4 + a2_25·x5 + a2_26·x6 = b2_2
a2_32·x2 + a2_33·x3 + a2_34·x4 + a2_35·x5 + a2_36·x6 = b2_3
...

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

k23 = −a2_32/a2_22

k24 = −a2_42/a2_22

k25 = −a2_52/a2_52

k26 = −a2_62/a2_22

И снова перепишем новые множители, теперь для уравнений с третьего по шестое:

a3_33 = a2_33 + k23·a2_23	a3_34 = a2_34 + k23·a2_24	a3_35 = a2_35 + k23·a2_25	...	b3_3 = b2_3 + k23·b2_2
a3_43 = a2_43 + k24·a2_23	a3_44 = a2_44 + k24·a2_24	a3_45 = a2_45 + k24·a2_25	...	b3_4 = b2_4 + k24·b2_2
...

И так далее, алгоритм вполне прозрачен. В итоге система уравнений принимает вид:

a11·x1 +	a12·x2 +	a13·x3 +	a14·x4 +	a15·x5 +	a16·x6	= b1
	a2_22·x2 +	a2_23·x3 +	a2_24·x4 +	a2_25·x5 +	a2_26·x6	= b2_2
		a3_33·x3 +	a3_34·x4 +	a3_35·x5 +	a3_36·x6	= b3_3
			a4_44·x4 +	a4_45·x5 +	a4_46·x6	= b4_4
				a5_55·x5 +	a5_56·x6	= b5_5
					a6_66·x6	= b6_6

Замечательная треугольная система, где все множители уже вычислены. Для ее решения остается пройтись по этой «лестнице» снизу вверх:

x6 = b6_6 / a6_66
x5 = (b5_5 − a5_56·x6) / a5_55
x4 = (b4_4 − a4_46·x6 − a4_45·x5) / a4_44
x3 = (b3_2 − a3_36·x6 − a3_35·x5 − a3_34·x4) / a3_33
x2 = (b2_2 − a2_26·x6 − a2_25·x5 − a2_24·x4 − a2_23·x3) / a2_22
x1 = (b1 − a16·x6 − a15·x5 − a14·x4 − a13·x3 − a12·x2) / a11

Становится примерно понятно, что методом Гаусса можно решить систему уравнений безумно большой размерности, но «на кончике пера» это получится уже с большим трудом. А вот с использованием современной вычислительной техники задача упрощается. Метод формализован, алгоритм решения получается несложный, а с объемом работы компьютер справится.

К слову, координаты максимума выпуклости колокола находятся так:
x0 = (2·b·d − c·f) / (c² − 4·a·b)
y0 = (2·a·f − c·d) / (c² − 4·a·b)

5. Оптимизация суммирования последовательности натуральных чисел

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

Но суммирование просто напрашивается на оптимизацию. При подготовке уравнения к решению множители при неизвестных вычисляются при помощи суммирования. Первое, что приходит на ум для программной реализации – цикл «for». И, пожалуй, трудно придумать ему замену при вычислении выражения, в котором участвует z, то есть, значения элемента массива. А вот те суммы, в которых есть только х и y, можно выразить несколько иначе. Ведь х и y – индексы массива, а значит, натуральные числа. А для вычисления суммы конечной последовательности натуральных чисел известны эквивалентные формулы:

1 + 2 + ... + n = n·(n + 1) / 2
1² + 2² + ... + n² = n·(n + 1)·(2·n + 1) / 6
1³ + 2³ + ... + n³ = n²·(n + 1)² / 4
1⁴ + 2⁴ + ... + n⁴ = n·(n + 1)·(2·n + 1)·(3·n² + 3·n − 1) / 30

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

В двумерном массиве x_i,j = i, y_i,j = j. В каждой строке значение х представлено рядом натуральных чисел от начального индекса до конечного. И в этом плане все строки массива одинаковы. Стало быть, суммы всех x_i, также, как и x_iⁿ, во всех строках одинаковы. Следовательно, сумма всех x_i,j в любой степени по всему массиву будет равна сумме x_i в соответствующей степени по всем элементам одной строки, умноженной на число строк.

То есть, для массива размерности m x n
∑(x) = ∑_j∑_i(x_i,j) = ∑_j(n·(n + 1) / 2) = m·n·(n + 1) / 2,
если индексация идет с единицы, и
∑(x) = (m + 1)·n·(n + 1) / 2,
если индексация идет с нуля. Соответственно, для y (при индексации с нуля) получается аналогично:
∑(y) = (n + 1)·m·(m + 1) / 2

Отдельно стоит обратить внимание на те суммы, где присутствуют и х, и y.
∑(x·y) = ∑_j∑_i(x_i,j·y_i,j)
Рассмотрим сначала вложенную сумму. К слову, сумма по i суммы по j равна сумме по j суммы по i, так что, в качестве вложенной будем рассматривать ту, которая в данный момент удобнее.
∑_i(x_i,j·y_i,j) = ∑_i(i·j)
Здесь j не зависит от i, а суммирование идет по i, стало быть, j можно вынести за знак суммы:
∑_i(i·j) = j·∑_i(i) = j·n·(n + 1) / 2
Получили выражение под внешней суммой:
∑_j(j·n·(n + 1) / 2)
Здесь n не зависит от j, и выражение с n можно вынести за знак суммы:
(n·(n + 1) / 2)·∑_j(j) = (n·(n + 1) / 2)·(m·(m + 1) / 2)

Аналогично и с другими выражениями, содержащими под общей суммой x и y, например:
∑(x³·y²) = (n²·(n + 1)² / 4)·(m·(m + 1)·(2·m + 1) / 6)

При вычислении суммы, содержащей z или ln(z) без цикла не обойтись. Но и здесь можно определить порядок суммирования по i или j, как удобно, вынести за знак вложенной суммы х или y в нужной степени, не зависящий от индекса суммирования, в результате операция умножения будет выполняться в m или n раз реже. Пример – ниже, в коде, в п. 7.

Следует учитывать, что формула вычисления суммы конечной последовательности натуральных чисел одинакова в случае, если ряд начинается с нуля или единицы, это очевидно. Если же требуется вычислить сумму последовательности натуральных чисел, начинающейся с числа k > 1 и заканчивающейся на n, необходимо из вычисленной суммы последовательности от 1 до n вычесть аналогичным образом вычисленную сумму последовательности от 1 до k. Здесь объем вычислений возрастает вдвое, и при программной реализации следует учесть число элементов последовательности n − k. При малых размерах организация вычисления через цикл может оказаться менее ресурсоемкой.

6. Замечания перед кодированием

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

При этом можно выполнить аппроксимацию массива, используя не весь массив, а лишь определенный его участок, скажем, квадратную область в центре или некоторую зону вблизи приблизительно определенных координат пика. Для этого необходимо всего лишь при вычислении сумм выполнять суммирование не от нуля до числа столбцов по x и от нуля до числа строк по y, а в определенных пределах, от x₀ − k до x₀ + k по x и от y₀ − p до y₀ + p по y.

При вычислении сумм используются значения максимальных (возможно, и минимальных) индексов массива. Так как они используются в том числе и в цикле «for» для суммирования z перебором массива, переменные для хранения этих индексов должны быть типа integer (скажем так, целого типа). Однако, при вычислении сумм x и y промежуточные значения могут выходить за пределы разрядной сетки даже и 64-разрядного целого. И неважно, что конечный результат будет присвоен переменной типа double, где он прекрасно поместился бы. Он будет сперва обрезан, а потом уже присвоен.

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

И еще. Каждый шаг алгоритма Гаусса похож на предыдущий. Здесь так и напрашивается рекурсия. Особенно при решении системы, в которой больше трех уравнений. Но если нет желания тратить время на качественную ее организацию, можно оправдать могучую технологию «Copy – Paste» тем, что рекурсия более ресурсоемка. Хотя, можно было бы написать красивую функцию для решения систем уравнений с возможностью задания числа уравнений в системе. Возможно, она даже уже кем-то написана...

7. Пример использования на Delphi

Вот пример использования методики – две процедуры на Delphi. Код не «причесан» и вряд ли может служить образцом хорошего программирования. Кроме того, для сокращения объема текста часть кода опущена в тех местах, где его не составит труда восстановить, используя материал теоретической части и имеющийся код. Это сделано во второй процедуре, и в местах пропусков вставлен комментарий {...и так далее...}.

7.1 Процедура аппроксимации плоскостью.

procedure Flat(Surf: pointer; out x1, x2, x3 : double); type // определяется тип-указатель на двумерный динамический массив. TArray = array of array of double; PArray = ^TArray; var // переопределим указатель на массив, чтоб узнать размерность Sur : PArray; // здесь будет размерность массива по х и y xm, ym : integer; // множители при неизвестных в системе уравнений // и правая часть системы, // индекс соответствует номеру уравнения a, b, c, d : array[1..3] of double; // накопитель для суммирования и коэффициенты для метода Гаусса buf, k12, k13, k23 : double; // счетчики цикла i, j : integer; begin Sur := PArray(Surf); // приведение типа указателя // выяснение размерности массива ym := Length(Sur^) - 1; xm := Length(Sur^[0]) - 1; // вычисление множителей при неизвестных в системе уравнений a[1] := (xm * (xm + 1) * (2 * xm + 1)) / 6 * (ym + 1); b[1] := xm * (xm + 1) * ym * (ym + 1) / 4; c[1] := xm * (xm + 1) / 2 * (ym + 1); a[2] := b[1]; b[2] := (ym * (ym + 1) * (2 * ym + 1)) / 6 * (xm + 1); c[2] := ym * (ym + 1) / 2 * (xm + 1); a[3] := c[1]; b[3] := c[2]; c[3] := xm * ym; // вычисление правой части системы уравнений d[1] := 0; buf := 0; for i := 0 to xm do begin for j := 0 to ym do buf := buf + Sur^[j, i]; d[1] := d[1] + buf * i; buf := 0; end; d[2] := 0; buf := 0; for i := 0 to ym do begin for j := 0 to xm do buf := buf + Sur^[i, j]; d[2] := d[2] + buf * i; buf := 0; end; d[3] := 0; for i := 0 to xm do for j := 0 to ym do d[3] := d[3] + Sur^[j, i]; // вычисление коэффициентов для умножения на 1-е уравнение // для преобразования 2-го и 3-го уравнений системы k12 := - a[2] / a[1]; k13 := - a[3] / a[1]; // вычисление коэффициента для умножения на 2-е уравнение // для преобразования 3-го уравнения системы k23 := - (b[3] + k13 * b[1]) / (b[2] + k12 * b[1]); // такой вариант вызывает ошибку переполнения типа integer. // исключительной ситуации не возникает, а результат неверный {x3 := ((d[3] + k13 * d[1]) + k23 * (d[2] + k12 * d[1])) / ((c[3] + k13 * c[1]) + k23 * (c[2] + k12 * c[1])); } // вычисление корней уравнения x3 := ((d[3] + k13 * d[1]) + k23 * (d[2] + k12 * d[1])); x3 := x3 / ((c[3] + k13 * c[1]) + k23 * (c[2] + k12 * c[1])); x2 := ((d[2] + k12 * d[1]) - (c[2] + k12 * c[1]) * x3) / (b[2] + k12 * b[1]); x1 := (d[1] - c[1] * x3 - b[1] * x2) / a[1]; end;

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

a) б)
Рисунок 2. Пример аппроксимации плоскостью

7.2 Процедура аппроксимации колоколом.

procedure Bell(Surf: pointer; out MaxX, MaxY : double; out x_1, x_2, x_3, x_4, x_5, x_6 : double); type TArray = array of array of double; PArray = ^TArray; var Sur : PArray; xm, ym : integer; // индексы множителей соответствуют номерам строки и столбца a, a2, a3, a4, a5, a6, k : array[1..6] of array[1..6] of double; // индексы правой части соответствуют номеру строки (уравнения) b, b2, b3, b4, b5, b6 : array[1..6] of double; buf : double; i, j : integer; // переменные для определения области аппроксимации // значения максимума, минимума и среднего значения массива max, min, half : double; // примерные координаты пика x0, y0 : integer; // индексы границ области аппроксимации x1_, y1_, x2_, y2_ : integer; // для цикла for x1, y1, x2, y2 : double; // для альтернативного вычисления t1 : double; begin // приведение типа указателя и выяснение размерности массива Sur := PArray(Surf); ym := Length(Sur^) - 1; xm := Length(Sur^[0]) - 1; // выяснение координат максимального элемента массива x0 := 0; y0 := 0; x1_ := 0; y1_ := 0; x2_ := 0; y2_ := 0; max := 0; for i := 0 to xm do for j := 0 to ym do if max < Sur^[j, i] then begin max := Sur^[j, i]; x0 := i; y0 := j; end; // вычисление расстояний от пика до середины склона // в четырех направлениях min := Sur^[y0, 1]; for i := x0 to xm do if min > Sur^[y0, i] then min := Sur^[y0, i]; half := max - ((max - min) * 0.5); for i := x0 to xm do if half > Sur^[y0, i] then begin x1_ := i; break; end; min := Sur^[y0, 1]; for i := x0 downto 0 do if min > Sur^[y0, i] then min := Sur^[y0, i]; half := max - ((max - min) * 0.5); for i := x0 downto 0 do if half > Sur^[y0, i] then begin x2_ := i; break; end; min := Sur^[1, x0]; for i := y0 to ym do if min > Sur^[i, x0] then min := Sur^[i, x0]; half := max - ((max - min) * 0.5); for i := y0 to ym do if half > Sur^[i, x0] then begin y1_ := i; break; end; min := Sur^[1, x0]; for i := y0 downto 0 do if min > Sur^[i, x0] then min := Sur^[i, x0]; half := max - ((max - min) * 0.5); for i := y0 downto 0 do if half > Sur^[i, x0] then begin y2_ := i; break; end; // переопределение границ области аппроксимации // - приведение к типу double x2_ := x2_ - 1; y2_ := y2_ - 1; x1 := x1_; y1 := y1_; x2 := x2_; y2 := y2_; // вычисление сумм - // это будут множители при неизвестных для исходного уравнения a[1, 1] := ((x1 * (x1 + 1) * (2 * x1 + 1) * (3 * Sqr(x1) + 3 * x1 - 1)) - (x2 * (x2 + 1) * (2 * x2 + 1) * (3 * Sqr(x2) + 3 * x2 - 1))) / 30 * (y1 - y2); a[1, 2] := (x1 * (x1 + 1) * (2 * x1 + 1) - x2 * (x2 + 1) * (2 * x2 + 1)) * (y1 * (y1 + 1) * (2 * y1 + 1) - y2 * (y2 + 1) * (2 * y2 + 1)) / 36; a[1, 3] := (Sqr(x1) * Sqr(x1 + 1) - Sqr(x2) * Sqr(x2 + 1)) * (y1 * (y1 + 1) - y2 * (y2 + 1)) / 8; a[1, 4] := ((Sqr(x1) * Sqr(x1 + 1)) - (Sqr(x2) * Sqr(x2 + 1))) / 4 * (y1 - y2); a[1, 5] := (x1 * (x1 + 1) * (2 * x1 + 1) - x2 * (x2 + 1) * (2 * x2 + 1)) * (y1 * (y1 + 1) - y2 * (y2 + 1)) / 12; a[1, 6] := (x1 * (x1 + 1) * (2 * x1 + 1) - x2 * (x2 + 1) * (2 * x2 + 1)) / 6 * (y1 - y2); a[2, 1] := a[1, 2]; a[2, 2] := ((y1 * (y1 + 1) * (2 * y1 + 1) * (3 * Sqr(y1) + 3 * y1 - 1)) - (y2 * (y2 + 1) * (2 * y2 + 1) * (3 * Sqr(y2) + 3 * y2 - 1))) / 30 * (x1 - x2); a[2, 3] := (Sqr(y1) * Sqr(y1 + 1) - Sqr(y2) * Sqr(y2 + 1)) * (x1 * (x1 + 1) - x2 * (x2 + 1)) / 8; a[2, 4] := (y1 * (y1 + 1) * (2 * y1 + 1) - y2 * (y2 + 1) * (2 * y2 + 1)) * (x1 * (x1 + 1) - x2 * (x2 + 1)) / 12; a[2, 5] := ((Sqr(y1) * Sqr(y1 + 1)) - (Sqr(y2) * Sqr(y2 + 1))) / 4 * (x1 - x2); a[2, 6] := (y1 * (y1 + 1) * (2 * y1 + 1) - y2 * (y2 + 1) * (2 * y2 + 1)) / 6 * (x1 - x2); a[3, 1] := a[1, 3]; a[3, 2] := a[2, 3]; { ... и так далее ... } // вычисление правой части системы // первое уравнение: b[1] := 0; buf := 0; for i := x2_ + 1 to x1_ do begin for j := y2_ + 1 to y1_ do buf := buf + Ln(Sur^[j, i]); b[1] := b[1] + buf * Sqr(i); buf := 0; end; { ... и так далее до b[6] ... } // шестое уравнение: b[6] := 0; for i := x2_ + 1 to x1_ do for j := y2_ + 1 to y1_ do b[6] := b[6] + Ln(Sur^[j, i]); // вычисление коэффициентов для умножения на 1-е уравнение // для преобразования 2, 3, 4, 5 и 6 уравнений k[1, 2] := - a[2, 1] / a[1, 1]; k[1, 3] := - a[3, 1] / a[1, 1]; k[1, 4] := - a[4, 1] / a[1, 1]; k[1, 5] := - a[5, 1] / a[1, 1]; k[1, 6] := - a[6, 1] / a[1, 1]; // вычисление новых коэффициентов для уравнений со второго и ниже // для второго уравнения: a2[2, 2] := a[2, 2] + k[1, 2] * a[1, 2]; a2[2, 3] := a[2, 3] + k[1, 2] * a[1, 3]; a2[2, 4] := a[2, 4] + k[1, 2] * a[1, 4]; a2[2, 5] := a[2, 5] + k[1, 2] * a[1, 5]; a2[2, 6] := a[2, 6] + k[1, 2] * a[1, 6]; // для третьего уравнения: a2[3, 2] := a[3, 2] + k[1, 3] * a[1, 2]; a2[3, 3] := a[3, 3] + k[1, 3] * a[1, 3]; { ... и так далее для остальных уравнений ... } // вычисление правой части для уравнений со второго и ниже b2[2] := b[2] + k[1, 2] * b[1]; b2[3] := b[3] + k[1, 3] * b[1]; b2[4] := b[4] + k[1, 4] * b[1]; b2[5] := b[5] + k[1, 5] * b[1]; b2[6] := b[6] + k[1, 6] * b[1]; // вычисление коэффициентов для умножения на 2-е уравнение // для преобразования 3, 4, 5 и 6 уравнений k[2, 3] := - a2[3, 2] / a2[2, 2]; k[2, 4] := - a2[4, 2] / a2[2, 2]; k[2, 5] := - a2[5, 2] / a2[2, 2]; k[2, 6] := - a2[6, 2] / a2[2, 2]; // вычисление новых коэффициентов для уравнений с третьего и ниже // для третьего уравнения: a3[3, 3] := a2[3, 3] + k[2, 3] * a2[2, 3]; a3[3, 4] := a2[3, 4] + k[2, 3] * a2[2, 4]; a3[3, 5] := a2[3, 5] + k[2, 3] * a2[2, 5]; a3[3, 6] := a2[3, 6] + k[2, 3] * a2[2, 6]; // для четвертого уравнения: a3[4, 3] := a2[4, 3] + k[2, 4] * a2[2, 3]; { ... и так далее для остальных уравнений ... } // вычисление правой части для уравнений с третьего и ниже b3[3] := b2[3] + k[2, 3] * b2[2]; b3[4] := b2[4] + k[2, 4] * b2[2]; b3[5] := b2[5] + k[2, 5] * b2[2]; b3[6] := b2[6] + k[2, 6] * b2[2]; { ... и так далее до преобразования шестого уравнения ... } k[5, 6] := - a5[6, 5] / a5[5, 5]; a6[6, 6] := a5[6, 6] + k[5, 6] * a5[5, 6]; b6[6] := b5[6] + k[5, 6] * b5[5]; // вычисление неизвестных x_6 := b6[6] / a6[6, 6]; x_5 := (b5[5] - a5[5, 6] * x_6) / a5[5, 5]; x_4 := (b4[4] - a4[4, 6] * x_6 - a4[4, 5] * x_5) / a4[4, 4]; x_3 := (b3[3] - a3[3, 6] * x_6 - a3[3, 5] * x_5 - a3[3, 4] * x_4) / a3[3, 3]; x_2 := (b2[2] - a2[2, 6] * x_6 - a2[2, 5] * x_5 - a2[2, 4] * x_4 - a2[2, 3] * x_3) / a2[2, 2]; x_1 := (b[1] - a[1, 6] * x_6 - a[1, 5] * x_5 - a[1, 4] * x_4 - a[1, 3] * x_3 - a[1, 2] * x_2) / a[1, 1]; // Вычисление координат максимума пика. // Если координаты (0, 0), значит, где-то ошибка. MaxX := 0; MaxY := 0; t1 := (Sqr(x_3) - 4 * x_1 * x_2); if t1 <> 0 then begin MaxY := (2 * x_1 * x_5 - x_3 * x_4) / t1; MaxX := (2 * x_2 * x_4 - x_3 * x_5) / t1; end; end;

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

a) б)
Рисунок 3. Пример аппроксимации колоколом

8. Используемая литература

1. Большая Советская Энциклопедия
2. Бронштейн И. Н., Семендяев К. А. Справочник по математике для инженеров и учащихся втузов. – 13-е изд., исправленное. – М.: Наука, Гл. ред.физ.-мат. лит., 1986. – 544 с.
3. Выгодский М. Я. Справочник по высшей математике – М.: Астрель, 2006. – 991,[1]с.:ил.