Унифицированный алгоритм дифференцирования и интегрирования дискретных числовых массивов Текст научной статьи по специальности «Математика»

Унифицированный алгоритм дифференцирования и интегрирования дискретных числовых массивов

В.В. Митюков, Ульяновское высшее авиационное училище гражданской авиации (институт),программист ОВТИ,

V. mityukov@gmail. сот

Введение

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

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

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

Обобщенный алгоритм

Предварительно принимается, что гладкое приближение заданных дискретных точек {х^ у^ на плоскости, осуществляется линейной аналитической моделью (например, степенным рядом):

/(х, Со ...С„) = Соф о(х) + С ф 1(х) + . . . + С„ф „(х) (1)

где С ^ - некоторые коэффициенты

ф¿(х) - базисные функции (линейно-независимые)

Условия интерполяции, то есть условия точного выполнения равенств /(х) = yi в точках хь (/' = 0, 1, .. т), приводят к системе линейных уравнений с квадратной матрицей плана Р (причем „ = т).

0 (х0 ) 1 ((Хо ) 0 ((Х1 ) 1((Х1 )

„(хо )\( Со Л = (Уо Л или Р с = у (2)

А)

V о

(Хт ) 1 (Хт )

1 ((Хт )) V Сп J

С1

У1

V ут J

При аппроксимации минимизируется сумма квадратов отклонений /(хк) от ук. на заданной системе точек хк , (к = 0, 1, .. т) (метод наименьших квадратов). Необходимое для этого условие экстремума приводит к системе нормальных уравнений с квадратной матрицей N размера п+1 х п+1 (причем, должно соблюдаться п < т):

N•0 = Ь (3)

Предлагается алгоритм [4], в котором вычисляемое значение / в точке х, т.е. / = / (х, Ск.), определяется из условия разрешимости расширенной однородной системы уравнений. Эта система получается из системы (2) или (3) переносом из ее правой в левую часть столбца у или Ь и добавлением к ней нижней строки:

Со ф о(х) + С ф 1(х) + . . + Сп ф п(х) -/= 0

(б

/о \( Со Л = (о Л

/1

/п

V о

(Х) 1(Х) - п(х) - /

С Сп

V! J

о

V о

или Н с = 0

(4)

где через О обозначена матрица Р или N через / значения -у,- или -Ь,.

Линейные преобразования Гаусса приведут квадратную матрицу Н размера г х г (п+2 х п+2) к треугольному виду и det(H) будет равен произведению ее диагональных элементов. Если в этом процессе не подвергать перестановкам нижнюю строку, то условие разрешимости [1] однородной системы (4), примет вид:

det(H) = det(Q)•(- / + к г г) = 0 (5) где к г г - линейная комбинация значений /, добавленная в нижний

о

правый элемент матрицы Н после преобразований Гаусса.

Поскольку заведомо det(Q) ф 0, что необходимо для существования единственного решения системы ( 2 ) или ( 3 ), то должно выполняться соотношение (- / + к г г) = 0 или к г г = /. Таким образом, алгоритм вычисления / состоит в обнулении правого нижнего элемента матрицы Н и последующего прямого хода алгоритма Гаусса без перестановок нижней строки. Накопленная в правом нижнем элементе матрицы Н линейная комбинация значений /¡, даст искомое значение /.

В случаях, когда серия вычислений производится при неизменных исходных дискретных данных (точки {х^ у}), алгоритм Гаусса можно не повторять для всей матрицы Н. Поскольку в этих случаях будет меняться только нижняя строка, то приведенные к треугольному виду строки матрицы Н, кроме нижней, следует сохранять в памяти и применять алгоритм Гаусса только к нижней строке.

Вместо набора „+1 коэффициентов С) ряда (1), в обобщенном алгоритме используется треугольная дискретная структура размера (п+2)-(п+3)/2. Такая неэффективная на первый взгляд реструктуризация вычисляемых данных, предоставляет следующие возможности:

• получать результаты, в том числе производные и интегралы, путем

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

• уменьшать потери точности в вычислениях, поскольку выполняется

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

Последующее обобщение

Универсальность обобщенного алгоритма позволяет распространять его применение на новые задачи. Можно видеть, что полученный алгоритм остается справедливым не только при вычислении в точке х значения /, но и при вычислении результата любой операции над рядом (1), сохраняющей его линейность относительно коэффициентов Ck (к = 0, 1,

.. „). В частности, к таким операторам относятся производная 1 в точке х, или интеграл ¥ на интервале [х-1 , х]. Для их вычисления достаточно заменить нижнюю строку системы (4):

Со ф о(х) + С ф 1 (х) + . . + С„ф п(х) -/= 0

на продифференцированную (6), или проинтегрированную (7) строку:

Г Г Г ¿"Г

Со ф о (х) + С1 ф 1 (х) + . . . + С„ф п (х) - 1 = 0 (6)

Со

/ф0( х) <х

+ Сп

(х) <х _

^ = 0 (7)

В принципе, если возникнет такая практическая необходимость, ничего не мешает повторять дифференцирование и интегрирование и включать в нижнюю строку системы (4) соотношения для производных различных порядков (р = 1, 2, 3...) и/или для повторных интегралов различной кратности (р = -1,-2,-3...).

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

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

п 1 ■

0

1 о о

V

1 - г2

2 0

1 - 1 - ¿3

2 3 0

1 - О2

2 0

0

1 о

1

3 1

3

О2

1 -1-*4 1 -1 - О

1 4 '

1

5 -

5Л

1 - о3

2- О

2-1

1 - ¿4 4 0

О

з - о2

3- 2- О

(8)

У

Чтобы получить первые 4 члена полинома Чебышева, а также результаты их двойного дифференцирования и интегрирования, нужно матрицу (8) умножить на матрицу (9) из чисел, вычисляемых по рекуррентной формуле

10 -10 1 о - 3 2 0 4

+

х -1

х -1

В случае полиномов Бернштейна, все элементы числовой матрицы уже зависят от ее порядка (количества членов) хотя алгоритм их вычисления не сложен. Например, для 2-х, 3-х и 4-х членов, матрицы будут следующими:

г

О

-1

V

1

-2 1

О 2 -2

О О

1/

< 1 О О О

—3 3 О О

3 —6 3 О

V—1 3 —3 1

1

1

Практическое опробование изложенных алгоритмов осуществлялось в офисной программе "MS Excel". Вводился набор дискретных данных, выбиралась система базисных функций, назначался метод приближения, задавались нужные значения констант m и п. Построение расширенной матрицы H, ее обработка и вычисление требуемых результатов производились программно. Для факторизации (LU-разложения) [2], [3] матрицы H использовался алгоритм Гаусса в модификации Холесского, с частичным выбором ведущего элемента. Необходимые для этого макросы и подпрограммы были реализованы на встроенном языке программирования VBA.

Литература

1. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Численные методы: Учебное пособие. - М.: Наука. Гл. ред. физ.-мат. лит., 1987. -600 с.

2. Форсайт Дж., Малькольм М., Моулер К. Машинные методы математических вычислений. Пер. с англ. - М.: Мир, 1980. -280 с.

3. Райс Дж. Матричные вычисления и математическое обеспечение. Пер. с англ. - М.: Мир, 1984. -264 с.

4. Митюков В.В. Обобщенный алгоритм и дискретная унифицированная структура для вычислительных задач. «Современные информационные технологии и ^-образование». Сборник докладов научно-практической конференции: учебно-методическое пособие. Под ред. проф. В.А. Сухомлина. -М.: ИНТУИТ.РУ, 2009. - с. 675-681