Перейти к содержимому

Интерполяционный полином Лагранжа в Excel

В рамках подготовки курса для бакалавров МФТИ я понял, что в моем блоге не так много заметок по использованию Excel в математике и физике. Каково же было мое удивление, когда я обнаружил, что книг по этой теме на русском языке буквально единицы. Ранее я опубликовал Вильям Дж. Орвис. Excel для ученых, инженеров и студентов. В заметке представлены три варианта нахождения интерполяционного полинома Лагранжа: таблица на листе Excel, функция VBA, функция листа Excel на основе REDUCE и LAMBDA. В заметке использованы материалы книги Алексея Васильева Числовые расчеты в Excel. Бумажная и электронная версии книги доступны на сайте издательства.

Рис. 1. Вычисление интерполяционного полинома по методу Лагранжа; чтобы увеличить изображение кликните на нем правой кнопкой мыши и выберите Открыть картинку в новой вкладке

Скачать заметку в формате Word или pdf, примеры в архиве (внутри файл Excel с поддержкой макросов)

Алгоритм

Задача интерполирования обычно решается для того, чтобы «восстановить» по набору дискретных данных аналитическую функциональную зависимость. Нередко в качестве функции, на основе которой строится интерполяционная зависимость, выбирают полиномиальные выражения. Предположим, имеется набор узловых точек х1, х2, …, хn и набор значений y1, y2, …, yn неизвестной функции в этих точках. Задача – построить зависимость f(x) такую, чтобы соответствующая кривая проходила через все точки (xk, yk) (k = 1, 2, …, n), т.е. необходимо, чтобы для всех k выполнялись соотношения f(xk) = yk.

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

Схема Лагранжа предполагает, что соответствующее полиномиальное выражение, построенное по точкам (х0, y0), (х1, y1), …, (хn, yn), ищется в виде

Полином имеет степень n, поскольку строится по n + 1 точке: индексация точек (хm, уm) начинается с нуля, а последний индекс равен n. Функции ϕm(х) являются полиномами степени n, причем такими, что в узловых точках xk имеют место соотношения: ϕmk) = 0, если k ≠ m, и ϕmk) = 1, если k = m. Эти полиномы вычисляются в виде произведений

для всех m = 0, 1, 2, …, n.

Метод Лагранжа в таблице на листе Excel

Реализуем метод Лагранжа в Excel для табулированной по 11 равноудаленным узловым точкам (на интервале значений аргумента от -2π до 2 π) функции

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

Ячейки D4:N4 содержат значения узловых точек, на основе которых мы создаем интерполяционный полином. Значения полинома в узловых точках отображаются в ячейках D5:N5. В ячейках D3:N3 указаны индексы узловых точек.

Основные вычисления выполняются в ячейках D7:N107 и, как результат, в ячейках С7:С107 вычисляются значения интерполяционного полинома в точках, которые указаны в ячейках А7:А107. В ячейках D7:N107 содержатся значения функций ϕm(х) для разных аргументов интерполяционного полинома х и разных индексов узловых точек m. В каждой строке диапазона D7:N107 находятся значения функций ϕm(х) для одного и того же аргумента х, но разных индексов m. В столбцах диапазона D7:N107 содержатся значения для разных аргументов х, и одного и того же индекса m. В диапазоне D7:N107 три типа формул: для левого D7: D107 и правого N7:N107 краев и остальных столбцов Е7:М107. С формулами можно ознакомиться в приложенном Excel-файле.

Выделите диапазон А7:С107 и постройте график. Чтобы добавить узловые точки, скопируйте диапазон D4:N5 в буфер обмена, выделите диаграмму, пройдите Главная –> Вставить –> Специальная вставка. Настройте параметры в окне Специальная вставка. Нажмите Ok. Отформатируйте вставленный ряд: отмените линию и добавьте встроенный маркер.

Рис. 2. График функции f(x) = sin(x)/(1 + x2), интерполяционный полином Лагранжа и узловые точки

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

Функция VBA

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

Пояснение кода

Функция ЛАГРАНЖ() имеет три аргумента: диапазон ячеек со значениями узловых точек (Рх), диапазон ячеек со значениями интерполируемой функции в узловых точках (Ру), а также аргумент, для которого вычисляется значение интерполяционного полинома Лагранжа (z). Функция возвращает числовой результат типа Double.

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

Алгоритм использует количество узловых точек. Явно этот параметр в аргументе функции ЛАГРАНЖ() отсутствует. Но он вычисляется, как количество ячеек в диапазоне, переданном первым аргументом функции ЛАГРАНЖ(). Поэтому мы объявляем переменную n = Рх.Count. Свойство Count для диапазона возвращает количество ячеек.

Основные вычисления производятся в блоке из вложенных условных операторов. Для использования в этих операторах объявляются две целочисленные переменные i и j. В переменной L накапливается полиномиальная сумма, а в переменной phi произведение в соответствии с формулой (2).

Перед началом выполнения вложенных операторов цикла переменной L присваивается значение 0. Индексная переменная i во внешнем цикле пробегает значения от 1 до n. В начале каждой итерации переменной phi присваивается значение 1. После этого последовательно запускаются два идентичных цикла. Основное различие между ними – диапазон изменения индексной переменной j. Для первого цикла она изменяется от 1 до i-1, а для второго цикла — от i+1 до n. Таким образом, при фиксированном значении i переменная j пробегает все значения от 1 до n, за исключением значения i. За каждую такую итерацию переменная phi сначала умножается на величину (x-Px.Cells(j).Value), а затем делится на величину (Px.Cells(i).Value-Px.Cells(j).Value). Здесь следует учесть, что Рх.Сеlls(индекс).Value – это значение ячейки с указанным индексом в диапазоне Рх (т.е. это значение узловой точки).

После того как значение переменной phi (для данного значения i) вычислено, командой L = L+phi * Ру.Cells(i).Value к полиномиальной сумме добавляем очередное слагаемое. Здесь Py.Cells(i).Value – ссылка на значение ячейки в диапазоне Ру со значениями табулированной функции. В итоге значение переменной L возвращается как результат функции (команда ЛАГРАНЖ = L).

Работа функции ЛАГРАНЖ()

Функцию ЛАГРАНЖ() можно использовать на рабочем листе. Рассмотрим новый пример. Ячейки А4:В9 содержат данные об узловых точках и значениях функции f(x) = x*exp(-x). При этом в ячейках А4:А9 указаны несколько неравномерно распределенных на интервале от 0 до 7 точек.

Рис. 3. Пользовательская функции для вычисления интерполяционного полинома Лагранжа

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

Вычисление полинома Лагранжа на основе функций REDUCE и LAMBDA

В декабре 2020 года Microsoft анонсировал функцию LAMBDA, которая позволяет определять пользовательские функции, написанные на языке формул Excel. А в июле 2021 г. объявил о создании новых функций, основанных на LAMBDA. Я недавно описал работу с LAMBDA и новыми функциями Excel. Не могу сказать, что написание сложных конструкций на основе этих функций проще, чем кода VBA, но как учебный пример, это весьма интересно.[1]

Код функции REDUCE

Описание работы функции REDUCE

Функция REDUCE работает, как обычная функция листа. Внутри ее могут располагаться ссылки на ячейки и динамические массивы. В отличие от функции LAMBDA, функция REDUCE не требует предварительного именования. Алгоритм работы функции REDUCE похож на алгоритм кода VBA в функции ЛАГРАНЖ(): один внешний цикл и один внутренний (в коде VBA два внутренних цикла). Внешний цикл перебирает все i значений диапазона В4:В9 (см. рис. 4), умножая на значения phi, определяемые для каждого i во внутреннем цикле. Во внутреннем цикле задается стартовое значение phi = 1. Далее для фиксированного i перебираются все j значений диапазона А4:А9, и для всех i ≠ j значение phi, полученное на предыдущем шаге, умножается на некое значение, а для i = j значение phi не изменяется. (Благодаря такой проверке, вместо двух внутренних циклов в коде VBA, здесь используется один.)

Посмотрим, как этот алгоритм реализован в коде функции REDUCE. Цель функции REDUCE – обработать массив, и вернуть одно число. Функции REDUCE имеет три аргумента. Первые два – начальное значение (0) и обрабатываемый массив ($B$4#).

Третий аргумент – функция обработки элементов массива

Рис. 4. Работа функции REDUCE

Функция LAMBDA принимает столько же параметров, сколько передает функция REDUCE:

L – накопитель; Py – массив. L – это то значение, которое вернется функцией после обработки всех элементов массива. Начальное значение L = 0. Это значение определено первым аргументом функции REDUCE. Ру – массив $B$4#. Функция LAMBDA накапливает значения L в элементе формулы L + REDUCE(… Эта запись аналогична более привычной для программистов L = L + REDUCE(…

Функция LAMBDA реализует внешний цикл. Она стартует со значения L = 0 и для всех элементов массива Ру (он же $B$4#) выполняет действие

Здесь реализована сумма произведений элементов Ру и рассчитанных коэффициентов, которые возвращаются внутренним циклом на основе REDUCE(…

Внутренний цикл также основан на REDUCE

Результат REDUCE – коэффициент для умножения на элемент массива Ру. Функция REDUCE принимает два аргумента

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

LAMBDA принимает два параметра от REDUCE: начальное значение phi = 1 и массив Рх = $A$4:$A$9. Функция LET не выполняет расчеты, а позволяет упростить восприятие этого фрагмента формулы (подробнее см. здесь). LET определяет три переменные – индекс элемента массива $A$4:$A$9, обрабатываемого на текущем шаге:

… индекс элемента массива $B$4#, переданного из внешнего цикла:

… и значение элемента массива $A$4:$A$9, соответствующего индексу Ind_y из внешнего цикла:

Расчет функции LET происходит внутри оператора ЕСЛИ:

Для всех элементов массива, где Ind_x ≠ Ind_y накапливается новое значение phi:

Как обычно, эта запись эквивалента phi = phi*(…)

Если же индексы равны Ind_x = Ind_y, оставляем значение phi без изменения.

Ссылка на ячейку листа (D4#) – это ссылка на динамический массив аргументов, для которых вычисляется полином Лагранжа. На рис. 4 аргументы размещены в столбце D.

Поскольку использование REDUCE и LAMBDA – это всё же не классический цикл, мне не удалось найти вариант, как обращаться к элементу массива по его индексу. Поэтому в переменных Ind_x и Ind_y реализовано обращение по значению элемента: ПОИСКПОЗ(Px;$A$4:$A$9;0). Такой подход накладывает ограничения на аргументы и значения функции в узловых точках: не должно быть повторений. Если хотя бы в двух узловых точках значения функций совпадают, функция ПОИСКПОЗ() отработает не корректно.

[1] Это оригинальный метод. Он не описан в книге Алексея Васильева «Числовые расчеты в Excel».

Добавить комментарий

Ваш адрес email не будет опубликован. Обязательные поля помечены *