39. Аппроксимация функции, заданной таблично, методом наименьших квадратов

Аппроксимация таблично заданных функций
методом наименьших квадратов

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

 

Формальная постановка задачи осуществляется следующим образом. Пусть на отрезке [a, b] задана одномерная сетка

hx = {xi  / xi = xi – 1 + hi, hi > 0, i = 1, 2, 3, …, n; x0 = a, xn = b},

в узлах xi которой заданы значения yi = f(xi), i = 0, 1, 2, …, n – соответствующие значения функции f(x).

Пусть также для аппроксимации табличных данных выбран некоторый класс функций j (х, c0, c1, c2, …, cm), m < n, где c0, c1, c2, …, cm – неопределенные коэффициенты, выбор значений которых позволяет определить конкретную функцию из выбранного класса. Требуется найти значения коэффициентов c0, c1, c2,…, cm, для которых выполнено условие

j (хi, c0, c1, c2, …, cm))2 Þ min.          (3.5.1)

 

Критерий (3.5.1), на основании которого осуществляется выбор значений коэффициентов c0, c1, c2, …, cm, является стандартным, используемым в МНК. Выбранные в соответствии с этим критерием значения коэффициентов позволяют определить среди множества функций конкретную функцию, наиболее согласованную с табличными (экспериментальными) данными или, иначе, обеспечивающую наилучшее среднеквадратическое приближение.

 

Как отмечалось выше, функция j (х, c0, c1, c2, …, cm) называется моделью, тогда коэффициенты c0, c1, c2,…, cm будем называть параметрами
модели или модельными. В дальнейшем ограничимся рассмотрением случая, когда модель линейно зависит от параметров и ее можно представить в виде

j (х, c0, c1, c2, …, cm) = c0j0(х) + c1j1(х) + c2j2(х) + … + cmjm(х). (3.5.2)

 

Модель вида (3.5.2) часто называют обобщенным полиномом. Здесь j0(х), j1(х), j2(х), jm(х) – множество так называемых базисных функций. Базисные функции могут быть как линейными, так и нелинейными функциями переменной x. Независимо от этого модель (3.5.2) остается линейной, поскольку она линейно зависит от модельных параметров c0, c1, c2, …, cm.


В качестве базисных могут быть выбраны, например, степенные функции

j0(х) = 1, j1(х) = х, j2(х) = х2, …, jm(х) = хm.

 

Тогда модель будет представлять собой полином степени m:

j (х, c0, c1, c2, …, cm) = c0 + c1х + c2х2 + … + cmхm.                          (3.5.3)

 

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

 

Итак, для линейной модели (3.5.2) требуется найти значения модельных параметров c0, c1, c2, …, cm, обеспечивающих выполнение условия (3.5.1). Эту задачу можно решить двумя способами.

 

Первый способ

С математической точки зрения поставленная задача является задачей безусловной минимизации функции нескольких переменных. Функция F(c0, c1, c2, …, cm) является квадратической одноэкстремальной, поэтому ее минимум можно искать, исходя из необходимых условий экстремума для функций нескольких переменных.

 

Запишем эти условия:

       (3.5.4)

Сократив на 2, преобразуем систему линейных уравнений (3.5.4) к виду

  (3.5.5)

Уравнения (3.5.5) называют нормальными. Матрица системы линейных уравнений (3.5.5) имеет вид

Г = (3.5.6)

Матрицу Г называют матрицей Грамма.(? так в методичке)

 

Введем обозначения:

ji = (ji(х0), ji(х1), ji(х2), …, ji(хn))т, i = 0, 1, 2, …, m,                   (3.5.7)

y = (y0, y1, y2, …, yn)т; c = (c0, c1, c2, …, cm)т.

Тогда матрицу Грамма можно записать в более традиционном виде

Г = ,                                    (3.5.8)

а систему (3.5.6) – в матричной форме

Гc = g.                                                                                                   (3.5.9)

 

Здесь (ji, jk), i, k = 0, 1, 2, …, mcкалярное произведение векторов ji и jk;

g = (g0, g1, g2, …, gm)т; gi = (y, ji), i = 0, 1, 2, …, m.                        (3.5.10)

Очевидно, что матрица Г симметричная. Можно показать, что она положительно определенная. Определитель матрицы Г (определитель Грамма) отличен от нуля в силу линейной независимости базисных функций j0(х), j1(х), j2(х), jm(х). Следовательно, система (3.5.4) имеет единственное решение, которое может быть найдено любым из известных методов решения систем с симметричной положительно определенной матрицей. В результате решения системы (3.5.4) будут найдены значения модельных параметров  и определена наилучшая в смысле МНК модель j (х,

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

 

Метод 1

Первый способ

1. На отрезке [a, b] задать одномерную сетку

hx = {xi / xi = xi – 1 + hi, hi > 0, i = 1, 2, 3, …, n; x0 = a, xn = b}

и значения yi = f(xi) в узлах сетки xi, i = 0, 1, 2, …, n.

2. Задать систему базисных функций j0(х), j1(х), j2(х), jm(х).

3. В соответствии с (3.5.7) рассчитать компоненты векторов ji,
i = 0, 1, 2, …, m.

4. Руководствуясь (3.5.8) и (3.5.10), вычислить элементы матрицы Г
и компоненты вектора
g.

5. Определить с = – вектор параметров модели (3.5.2) как решение системы линейных уравнений (3.5.9).

6. Модель j (х,  построенная согласно (3.5.2), является наилучшей в смысле МНК. Процесс завершен.

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

Второй способ

Поскольку модель (3.5.2)

j (х, c0, c1, c2, …, cm) = c0j0(х) + c1j1(х) + c2j2(х) + … + cmjm(х)

должна соответствовать yi = f(xi), i = 0, 1, 2, …, n – данным, представленным в узлах сетки

hx = {xi  / xi = xi – 1 + hi, hi > 0, i = 1, 2, 3, …, n; x0 = a, xn = b},

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

j (хi, c0, c1, c2,…, cm) = yi, i = 0, 1, 2, …, n.

Эта система уравнений является системой линейных уравнений, поскольку модель j (хi, c0, c1, c2, …, cm) линейна относительно модельных параметров c0, c1, c2,…, cm. Запишем ее в детальной форме:

                (3.5.11)

Матрица системы (3.5.11) – прямоугольная n×m матрица:

А =                                            (3.5.12)

Система линейных уравнений (3.5.11) является переопределенной
с матрицей
A, имеющей полный столбцовый ранг, так как по предположению m < n, а базисные функции j0(х), j1(х), j2(х), jm(х) линейно независимы. Как известно, для переопределенных систем линейных уравнений ищется так называемое обобщенное решение, которое обеспечивает минимум суммы квадратов невязок

что соответствует условию (3.5.1).

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

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

АтА c = Ат y.                                       

Поскольку в обозначениях (3.5.7) и (3.5.8) Атy = g, а АтА = Г, полученная система эквивалентна системе (3.5.9).

Так как первая трансформация Гаусса ухудшает обусловленность системы линейных уравнений, то применение такого подхода не всегда целесообразно. В ряде случаев систему (3.5.11) предпочтительнее решать с помощью любого другого из известных методов, например с использованием ортогональных преобразований, что обеспечит большую устойчивость вычислительного процесса. В последнем случае процесс построения наилучшей в смысле метода наименьших квадратов аппроксимации таблично заданной на отрезке [a, b] функции можно представить в виде следующего алгоритма.

 

Метод 2

1. На отрезке [a, b] задать одномерную сетку

hx = {xi / xi = xi – 1 + hi, hi > 0, i = 1, 2, 3, …, n; x0 = a, xn = b}

и значения yi = f(xi) в узлах сетки xi, i = 0, 1, 2, …, n.

2. Задать систему базисных функций j0(х), j1(х), j2(х), jm(х).

3. Вычислить jk(хi), k = 0, 1, 2, …, m; i = 0, 1, 2, …, n – значения базисных функций в узлах сетки hx и в соответствии с (3.5.12) сформировать матрицу A.

4. Определить с = (– вектор параметров модели (3.5.2) как обобщенное решение системы линейных уравнений (3.5.11). Решать систему (3.5.11) можно помощью одного из методов, использу-ющих ортогональные преобразования.

5. Модель j (х, , построенная согласно (3.5.2), является наилучшей в смысле МНК. Процесс завершен.

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

1. В случае когда m – число модельных параметров будет равно n – количеству узлов сетки hx, наилучшая в среднеквадратическом смысле модель j (х,  будет интерполяционной.

2. Метод 1 может быть предпочтительней метода 2 при m << n.
В этом случае преимущество первого метода будет заключаться в небольшой размерности системы линейных уравнений (3.5.9) по сравнению с системой (3.5.11).

3. Метод 2, очевидно, будет предпочтительней метода 1 при
не очень хорошей или плохой обусловленности матрицы (3.5.12). В этом случае вычислительная устойчивость процесса будет в большей степени гарантирована при решении системы (3.5.11) без применения первой трансформации Гаусса.

 

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

 

В целом следует отметить большую практическую значимость решения проблемы аппроксимации функций и МНК.