Как правильно дифференцировать дискретные функции (Часть 2. Все-таки, МКЭ?) / Хабр

В части первой обсуждались алгоритмы вычисления первой производной дискретной функции (функции, заданной массивом {аргумент: значение}, или массивом значений в узлах на равномерной сетке). Были протестированы функции Python numpy.diff и numpy.gradient (других в Python так и не нашлось), и был предложен алгоритм improved, который, как выяснилось, работает хорошо, но не всегда. При этом проблемы возникают на концах отрезка, на котором определена дифференцируемая функция. Во второй части поработаем с методом конечных элементов (МКЭ, Finite Elements Method), и попробуем построить алгоритм вычисления первой производной, наилучший по точности, и, по возможности, без “слабых мест”.

Во избежание разногласий (а они были в комментариях к части первой) еще разок об объекте исследования. В части первой было сделано одно не совсем корректное допущение. Все исследуемые массивы "визуально" были представлены с точностью до 3-х знаков после запятой. Однако в расчетах использовались точные значения, сгенерированные функциями numpy.exp(x) или x**3. Точность, естественно, была ограничена ошибками машинного округления. Поскольку в этой части повествования мы действуем более аккуратно, то рассматриваться будут дискретные функции, заданные именно массивами с точностью до 4-го знака после запятой. Например, таким:

$[1.0000, 1.1052, 1.2214, 1.3499, 1.4918, …, 2.7183]\quad\quad(1)$

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

Кроме массива (1) будут исследованы еще несколько других массивов. На них можно посмотреть здесь:

Массивы (дискретные функции) для исследования

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

Первая производная для оценки точности результатов:

Сгенерировано функцией на отрезке $[0, \pi/2]$ :

Первая производная для оценки точности результатов:

Про дискретную функцию (1) мы также знаем, что эта функция задана на отрезке , разделенном на части с равномерным шагом h=0.1. Алгоритм numpy.diff мы из рассмотрения исключим, как “наихудший” из-за несохранения размерности массива данных, и соревнование будет проводиться между остальными двумя ( numpy.gradient и improved), а также новыми алгоритмами, построенными на МКЭ.

Немного теории. Нудновато, но полезно

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

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

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

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

$<f, g> = [\int_{a}^{b} f(x)g(x) \,dx]\quad\quad(2)$

, которое определяет меру (модуль функции) :

$\parallel f \parallel^2 \quad = \quad<f, f> \quad= \int_{a}^{b} f(x)^2 \,dx$

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

$e_i(x), \quad i=0,\quad … , \quad n$

Изначальная задача - найти набор констант

$a_i, \quad i=0,\quad … , \quad n$

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

$\{a_i\} => min \quad \delta \quad = \quad min \int_{a}^{b} [f(x)- \sum \limits_{i=0}^n a_i*e_i]^2 \,dx$

Для этого необходимо равенство нулю всех частных производных

$\frac{\partial \delta}{\partial a_i}, \quad i=0,\quad... ,\quad n \quad\quad (3)$

Это и есть система n+1 линейных уравнений для нахождения n+1 неизвестных величин . В простейшей версии МКЭ для одномерного случая в качестве базисного набора функций выбираются “функции-крышки” - кусочно-линейные функции , такие, что

$e_i(x)=1 \quad при \quad x=x_i\quad и \quad e_i(x)=0 \quadпри\quad x\neq x_i$

В этом случае коэффициенты суть не что иное как приближенные значения функции в узлах сетки , что можно обозначить как $f(x_i) \quad \approx \quad a_i$ , а функция аппроксимируется конечно-линейной функцией - линейной комбинацией базисных функций:

$f(x) \quad \approx \quad\sum \limits_{i=0}^n a_i*e_i$

Важное замечание, в том числе, по терминологии

В дальнейшем изложении кусочно-линейная функция иногда будет называться вектором $\overline{a}$ . Этот вектор состоит из коэффициентов разложения функции по базису $\{e_i\}$ .

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

Опуская выкладки, итоговую систему линейных уравнений (3) можно представить в следующем виде:

$G * \overline{a} = \overline{<e, f>}\quad\quad (4)$

где G - квадратная матрица с компонентами $g_{ij}$ :

$g_{ij} = <e_i, e_j>$

где i, j = 0, …, n, вектор $\overline{a}$ - есть вектор-столбец из неизвестных величин , а вектор $\overline{<e, f>}$ -есть вектор-столбец из значений , где i = 0, …, n.

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

$\overline{a} = \overline{<e^*, f>}$

, где $\overline{<e^*, f>}$ - вектор-столбец из значений , и вектор $\overline{e^*}$ (эквивалентный набору векторов $\{e^*_i\}$ ) называется “сопряженным базисом” и находится по формуле:

$\overline{e^*} = G^{-1} * \overline{e}$

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

На практике матрицу G обычно явно не обращают, поскольку это очень затратная операция. Полезное свойство матрицы G в нашем (одномерном) случае - трехдиагональность, т.е. все ее ненулевые элементы расположены на главной диагонали (i = j) и на двух соседних диагоналях (| i - j | = 1). При этом система уравнений (4) решается прогонкой, с количеством операций O(n).

Функции-компоненты сопряженного базиса также обычно не вычисляют, однако ниже мы впрямую встретимся с ними, что лично мне кажется достаточно эстетичным.

МКЭ способ 1 ("лобовой")

Итак, начинаем применять теорию на практике. В прошлом дискретная, а теперь кусочно-линейная функция задана массивом значений (1) на отрезке с шагом h = 0.1, т.е. набором значений . При этом:

$f(x)= \sum \limits_{i=0}^n f_i*e_i \quad\quad(5)$

где - элементы массива (1), а - "функции-крышки". Необходимо найти ее первую производную, или, если точнее, наилучшую аппроксимацию первой производной, “натянутую” на базис $\{e_i\}$ .

Важный комментарий для понимания алгоритма

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

Обозначим первую производную функции как $f^{'}$ , и будем искать ее аппроксимацию на множестве линейных комбинаций функций-крышек, т е. в виде

$f^{'}(x) \quad \approx \quad\sum \limits_{i=0}^n a_i*e_i$

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

$\delta \quad = \quad \int_{a}^{b} [f^{'}(x)- \sum \limits_{i=0}^n a_i*e_i]^2 \,dx \quad\quad (6)$

Раскрываем квадратные скобки под интегралом, “выкидываем” первое слагаемое, поскольку оно не зависит от , и при дифференцировании все равно пропадет, затем ко второму слагаемому применяем интегрирование по частям, после этого система уравнений (3) (про равенство нулю всех частных производных $\delta$ по принимает вид:

$G * \overline{a} = f_n- f_0 -\overline{<{e^{'}}, f>} \quad\quad(7)$

где и - суть первый и последний элементы массива (1), соответственно, а $\overline {<e^{'}, f>}$ - вектор из скалярных произведений первых производных базисных функций $e^{'}_i = (e_i)^{'}$ и исходной функции . Поскольку базисные функции кусочно линейны, их производные $e^{'}_i$ будут кусочно-постоянны, и легко интегрируются.

Конкретный вид матрицы жесткости G и правой части уравнения (7) здесь:

$G = \begin{vmatrix} h/3 & h/6&0&... & 0& 0 & 0\\ h/6 & h/3 & h/6 &0&...& 0 & 0\\0 &h/6&h/3&h/6&0&...&0\\...&...&...&...&...&...&...\\...&...&0&h/6&h/3&h/6&0\\...&...&...&0&h/6&h/3&h/6\\...&...&...&0&0&h/6&h/3\end{vmatrix}$

Вектор-столбец в правой части уравнения (7):

$\begin{vmatrix} (f_1-f_0)/2\\ (f_2-f_0)/2\\...\\(f_{k+1}-f_{k-1})/2\\...\\(f_n-f_{n-2})/2\\(f_n-f_{n-1})/2\end{vmatrix}$

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

Рис.2. Сравнение алгоритмов на левом конце отрезка

Рис.3. Сравнение алгоритмов на правом конце отрезка

В середине отрезка у всех 3-х алгоритмов точность примерно одинакова, а проблемы, как обычно, на концах. При этом FEM устойчиво показывает результаты лучше, чем numpy.gradient, но хуже, чем improved.

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

Рис.4 Кубическая парабола на левом конце отрезка

На левом конце отрезка FEM практически точно ложится на тестовое решение. Вблизи правого конца отрезка "расклад сил" похож на Рис.3.

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

МКЭ способ 2 ("продвинутый")

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

$\delta \quad = \quad \int_{a}^{b} [f^{''}(x)- \sum \limits_{i=0}^n a_i*e_i]^2 \,dx \quad\quad (8)$

где $f^{''}(x)$ - вторая производная функции , которая, в свою очередь, представлена выражением (5). Казалось бы, вторая производная кусочно-линейной функции - тождественный (кроме точек разрыва) ноль. Ан нет, все не так просто. Повторяем процедуру так же, как в прошлом разделе. Раскрываем квадратные скобки под интегралом, “выкидываем” первое слагаемое, поскольку оно не зависит от , и при дифференцировании все равно пропадет, ко второму слагаемому применяем интегрирование по частям, после этого система уравнений (3) (про равенство нулю всех частных производных функции ошибки $\delta$ по переменным принимает вид (аналогично уравнению (7)):

$G * \overline{a} =f^{'}_n - f^{'}_0 - \overline{<e^{'}, f^{'}>} \quad\quad(9)$

где $f^{'}$ - первая производная функции , т.е. кусочно-постоянная функция, которая легко интегрируется в произведении с производными базисных функций $e^{'}_i$ - тоже кусочно-постоянными функциями, а $f^{'}_0$ и $f^{'}_n$ - значения первой производной на левом и правом концах отрезка , соответственно. При всем этом -

значения $f^{'}_0$ и $f^{'}_n$ заранее неизвестны и пока непонятно, откуда их взять.

Какая при этом получается система линейных уравнений? Матрица жесткости G остается такой же, как в предыдущем разделе (см. спрятанный текст выше), а вектор в правой части уравнения (9) (обозначим его $\overline{b}$ ) заслуживает подробного рассмотрения. Его можно представить в виде трех слагаемых:

$\overline{b} = \overline{b_{base}} -f^{'}_0*\overline{b_0}+f^{'}_n*\overline{b_n}$

Подробный вид слагаемых-векторов можно посмотреть здесь:

$\overline{b}=\begin{vmatrix} (f_1-f_0)/h - f^{'}_0\\ (f_2-2*f_1+f_0)/h\\...\\(f_{k+1}-2*f_k + f_{k-1})/h\\...\\(f_n-2*f_{n-1}+f_{n-2})/h\\f^{'}_n - (f_n-f_{n-1})/h\end{vmatrix}$ $\overline{b_{base}} = \begin{vmatrix} (f_1-f_0)/h \\(f_2-2*f_1+f_0)/h\\...\\(f_{k+1}-2*f_k + f_{k-1})/h\\...\\(f_n-2*f_{n-1}+f_{n-2})/h\\ - (f_n-f_{n-1})/h\end{vmatrix}$ $\overline{b_0} = \begin{vmatrix} 1\\0\\...\\0\\...\\0\\0\end{vmatrix} \quad\quad \overline{b_n}=\begin{vmatrix} 0\\0\\...\\0\\...\\0\\1\end{vmatrix}$

А значит, решение системы системы уравнений (9), представленной в виде:

$G*\overline{a} = \overline{b} \quad\quad(9.1)$

в свою очередь, представимо в виде:

$\overline{a} = \overline{a_{base}} -f^{'}_0*\overline{a_0}+f^{'}_n*\overline{a_n} \quad\quad(10)$

где $\overline{a_{base}}, \quad\overline{a_0}, \quad\overline{a_n}$ - решения 3-х систем линейных уравнений с одной и той же матрицей жесткости G и векторами в правой части $\overline{b_{base}}, \quad\overline{b_0}, \quad\overline{b_n}$ , соответственно. Эти три вектора-решения являются представлениями кусочно-линейных функций, “разложенных” по базису $\{e_i\}$ , причем компоненты каждого вектора являются коэффициентами при базисных "функциях-крышках" в этом разложении. Вектора ( т.е. функции) $\overline{a_{base}}, \quad\overline{a_0}, \quad\overline{a_n}$ для нашего исследуемого массива (1) представлены на рис.5.

Рис.5. Решение "base" и две функции из сопряженного базиса

Зеленый и красный графики почти сливаются на левом конце отрезка, поскольку разница между ними порядка 1, а масштаб рисунка по вертикали значительно больше (значение в нуле порядка 40). Заметим, что вектора $\overline{a_0}$ и $\overline{a_n}$ - не что иное как представление двух функций из сопряженного к $\{e_i\}$ базиса, а именно, функций и . Вот мы и встретились с сопряженным базисом, как было обещано в начале публикации.

Задача нахождения вектора $\overline{a}$ из выражения (10) поначалу выглядит достаточно сложной. Необходимо провести оптимизацию по двум параметрам $f^{'}_0$ и $f^{'}_n$ , причем цель оптимизации (типа минимизации какой-то новой функции ошибки) пока даже не сформулирована. Но если внимательно посмотреть на Рис.5, то можно заметить, что все три вектора из (10) имеют ярко выраженной "пилообразный" характер. Поэтому сформулируем:

Цель оптимизации - вычисление вектора $\overline{a}$ , соответствующего "наименее пилообразной" функции, т.е. наиболее гладкой аппроксимации $f^{''}$

На рис.5 хорошо видно, что для обоих участвующих в процессе сопряженных вектора $\overline{a_0}$ и $\overline{a_n}$ вся их пилообразность “сосредоточена” вблизи концов отрезка - вблизи левого конца отрезка для вектора (функции) с номером 0 (зеленый цвет) и вблизи правого конца отрезка для вектора (функции) с номером n (синий цвет). Сформулируем две гипотезы:

Гипотеза 1. Оптимизацию выражения (10) по двум параметрам $f^{'}_0$ и $f^{'}_n$ можно проводить независимо по каждому параметру отдельно. Первый будет отвечать за “гладкость” искомого вектора $\overline{a}$ ( т.е. функции $f^{''}$ ) слева, а второй - справа.

Гипотеза 2. В качестве “меры пилообразности” можно взять угол между первым и вторым звеньями ломаной $\overline{a}$ для нахождения оптимального значение параметра $f^{'}_0$ и угол между последним и предпоследним звеньями ломаной $\overline{a}$ для нахождения оптимального значения параметра $f^{'}_n$ . Вычислять будем косинус угла между звеньями, который, естественно, должен быть максимален в обоих случаях. Оптимизацию (поиск максимума косинуса) можно проводить разными хорошо известными методами, мы на этом останавливаться не будем.

Для оконгчательного решения поставленной задачи осталось по узловым значениям второй производной ( а это - вектор, он же набор узловых значений $\overline{a}$ ) восстановить набор узловых значений первой производной, например, так:

$f^{'}_{i+1} = f^{'}_i+h*(a_{i+1}+a_i)/2,$

где $f^{'}_0$ определено как значение первого параметра оптимизации. Оказалось, что эти две гипотезы прекрасно работают. Посмотрим на результат (solution) и сравним его с альтернативами. Сначала для экспоненты (ой, обещал забыть название функции):

Рис.6. Решение Solution соревнуется только с Improved и побеждает

Теперь для кубической параболы (ой, опять вспомнил название):

Рис.8. На левом крае Solution практически совпадает с FEM, оба побеждают

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

По ��екрету сообщаю, что сравнил результаты для функции $f=sin(x), \quad x\in[0,\pi/2]$ .

Результат здесь:

Слева: Solution опять соревнуется только с FEM

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

Наверное, в третью часть также добавлю Python-листинг. Пока код существует в Jupiter Notebook и требует доведения до презентабельного вида. Ну и кое-какие мелочи еще требуется доработать.

Спасибо за внимание. Буду благодарен за комменты и отзывы.