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

Митюков В.В.

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

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

КЛЮЧЕВЫЕ СЛОВА:

Метод гладкого приближения, обобщение, дискретное множество,

числовой ряд.

АНОТАЦИЯ:

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

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

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

Уравнения

Представляемый алгоритм гладкого приближения заданных точек {хг-у,} основывается на линейной модели, составленной из аналитически вычисляемых фрагментов - базисных функций (например, из членов степенного ряда или из членов ряда Чебышева, Фурье):

у(х) = Со • Ф0 (X) + С • %1 (х) +... + Сп •%„ (X) ( 1 )

где С j - искомые коэффициенты (/ = 0, 1, .. п);

Ф/(х) - базисные функции (аналитические и линейно-независимые).

Результаты операций дифференцирования и интегрирования линейной модели (1), остаются линейными относительно коэффициентов С/ , в виде:

у'(х) = Со • %0 (х) + С • Ф (х)+... + Сп • %п (х) ( 2 )

X

X

X

Y(х_ь х) = С0 • :%0(х) • dx + С1 • (х) • ёх +... + Сп • /фп(х) • ёх

х_

х _

х _

(3)

здесь [х-1, х] - интервал интегрирования.

Следует напомнить, что в вычислительной практике неизвестные коэффициенты С) в основном определяются из условий интерполяции или из условий аппроксимации [1], [2].

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

ЧоК) Ф1СО " Фп(x0)" 9 Со 6 9 уо 1

Фо(X1) %1(X1) " Фп(X1) "1 = у.

чФо(^ ) %1(^ ) " Фп(^ ), С п 8 ут 5

или Р с = у

( 4 )

Если имеются измеренные значения наклонов касательных уг- и/или значения интегральных площадей У, на некоторых подынтервалах (например [х^ Xi ]), то нет препятствий для пополнения системы уравнений (4), строками (равенствами) у'(х;) = у' и У(Xk, ) = У;, соответствующими (2), (3) [3]. Тогда линейную систему (4) с такими добавленными строками, следует записать в другом, более общем виде:

р00 р01 " ' р0п р10 р11 " ' р1 п

9 со 6 9 г 6 ^ о

С1 = 71

С 8^ п 7 8 т 5

или

Р с = z

( 5 )

8 рт 0 рт 1 " ' р

представляют собой значения базисных функций $ j (х,),

Здесь р ,

(или производных ф(х;), или интегралов /ф j (х) • ёх в добавленных

X ;

строках) .

Через z, обозначены значения у, (а также у', или У, ).

Размерность т системы (5), разумеется, будет превышать число точек х, на число дополнительных строк.

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

При аппроксимации минимизируется сумма квадратов отклонений у(х,) от у,. на заданной системе точек х, , (, = 0, 2, .. т) (метод наименьших квадратов [1]). Такой критерий позволяет использовать аналитический подход:

Для получения квадрата длины вектора невязки системы (5)

Аz = Р • с - z его нужно умножить скалярно самого на себя:

|| Аъ||2 = Аът • Аъ = (Рс - ъ)т • (Рс - ъ) Учитывая, что транспонирование сумм и произведений матриц А и В подчиняется следующим законам:

(А + В)т = Ат + Вт и (А-В)т = ВтАт

можно продолжить более подробно: (Рс - ъ)т<Рс - г) = ст Рт Рс

ст Р

ът Рс

т

Необходимым условием минимума скалярной функции а ъ| 2 является равенство нулю всех ее частных производных по векторному аргументу с :

& II А г ||

& с

= 0

или

Рт Рс + ст Рт Р - Рт 2 - ът Р = 0 .

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

21Рт Рс -21Ртъ

Это приводит к системе квадратной матрицей размера развернутый вид :

= 0 или линейных п+1 х п+1,

Рт Р с = Рт ъ

нормальных уравнений с принимающей следующий

Роо Рю Рог Рп

* * * рт 0

* " рт 1

Ро п Р

1 п

Р

Р00 Р01 Р10 Р11

Р 0 п

• Р1 п

тп 5 8

Рт 0 Р

т 1

Р

С

С1

с

Р00 Р10 Р01 Р11

*" Рт 0 * Рт 1

тп5 8сп 5 8Р0п Р1 п

Р

70 7Л

7

т п5 8 т 5

( 6 )

или в матрично-векторной записи N с = Ь

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

Алгоритм

В задачах аппроксимации прямого решения систем линейных уравнений (5) и (6) стараются избегать из-за того, что они часто получаются плохо обусловленными [1], [2]. Поэтому существуют другие способы, основанные на специально сконструированных базисных функциях [1]. Например, при интерполировании полиномами, путем разбивки и перегруппировки членов степенного ряда можно получить такие базисные функции как полиномы Лагранжа, в которых коэффициенты Си принимают значения равные заданным значениям уи Другими способами интерполяции являются полиномы Ньютона, итерационный процесс Эйткена и т.д..

В данном случае, для достижения универсальности используется

другой способ получения результата (произвольного значения у(х), производной у'(х), интеграла У(хх) ), далее обозначаемого через f Как показано в [3], сначала полученную систему уравнений (5) или (6) требуется расширить до однородного вида. Затем применяется единый подход к вычислению результатов, вытекающий из условия существования нетривиального решения однородной системы с матрицей Н, то есть из условия det (Н) = 0. Однородная система получается из системы (5) или (6) с квадратной матрицей Р или N переносом из ее правой в левую часть столбца z или Ь и добавлением к ней нижней строки, соответствующей (1), (2) или (3), общий вид которой :

со■ Ро(х) + с1-Р1(х) + ••• + сп■ Рп(х) -/ = 0 Таким образом:

-У-ь»)

% У-, -У-ь,

№5 - ■рп(х) -/ J

1

^ о ^ о

о о

нпн Н с = О

или Н с = 0 ( 7 )

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

det(H) = det(p) • (-/ + h г) = 0 или det(H) = det(N) • (-/ + h ь) = 0 где h 7 или h ь - линейная комбинация значений -2 i или -Ь i , добавленная в нижний правый элемент матрицы Н после преобразований Гаусса.

Поскольку заведомо det(P) \ 0 (или det(N) \ 0), что необходимо для существования единственного решения системы (5) или (6), то должны выполняться соотношения (-f + h 2) = 0 или (- f + h ь) =0, означающие, что / = h 2 или f = h ь .

Таким образом, алгоритм вычисления f состоит в обнулении правого нижнего элемента матрицы Н и последующего прямого хода алгоритма Гаусса без перестановок нижней строки. Накопленная в правом нижнем элементе матрицы Н линейная комбинация значений -2 или -Ьь определит искомое значение результата f .

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

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

виде [4],

¿0 -г,

— 7 \ Л/

-10 0 0 0-100 0 0-10 0 0 0 -1

(7 \ ¿0

Z

\ п/

11ГП1

■Ьо -¿1

-Ь

П/

f-\ о о 0-100 0 0-10 ч0 О 0 -1,

о Ï (ъ^

¿4

\bnj

то после преобразований Гаусса, результат будет представлен в виде линейной комбинации исходных значений zt или b i :

Wo(x)• Zo + Wi(x)• Zi + ... + wn(x)• zn или Wo(x)• b0 + Wi(x)• b + ... + wn(x)-bn В частности, таким путем можно генерировать различные квадратурные формулы (их веса и узлы), как традиционные (Ньютона-Котеса, Гаусса) [4], так и совершенно новые, для произвольного расположения узлов и выбора базисных функций.

Если же в развернутом виде представить нижнюю строку матрицы Н,

■ 1 0 о ^

(р0(х) р,{х) р„(х)'=(р0(х) рл(*)|-' 0 1

о о ■■■ 1 ^

то после преобразований Гаусса, результат будет представлен в виде линейной комбинации базисных функций (ф;(х), или % (х), или же

j (x)-dx ) , т.е. : Q- Ро(х) + Cr л(х) + •••+ Cn - pn (x) .

x _

Наиболее же общий, развернутый вид результата будет выглядеть следующим образом:

\р0(х) р^Х) ■■■ рп(г)).

f/i d ---d ^ a 00 ы 01 "on

¿Î j Q Yi " " " ¿Î

In

dn g dfi 1 ■ ■ ■ dnn

-zj-b, -ZnhK

Эскизное опробование изложенного алгоритма проводилось в офисной программе "MS Excel". Вводился набор дискретных данных, выбиралась система базисных функций, назначался метод приближения, задавались нужные значения констант тип. Построение расширенной матрицы H, ее обработка и вычисление требуемых результатов производились программно. Необходимые для этого макросы и подпрограммы реализованы на встроенном языке программирования VBA. Для полученных результатов предусмотрен вывод в таблицы и отображение на графиках.

Выводы

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

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

интегральных).

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

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

Литература

1. Каханер Д., Моулер К., Нэш С. Численные методы и математическое обеспечение. Пер. с англ. - М.: Мир, 1998. -575 с.

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

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

4. Митюков В.В. Задача унификации квадратурных формул численного интегрирования. «Современные информационные технологии и ИТ-образование». [Электронный ресурс] / Сборник научных трудов VI Международной научно-практической конференции / Под ред. В.А. Сухомлина. - Москва: МГУ, 2011. - Т. 1. - с. 546-551