Среднеквадратичное приближение функций

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

Теория


Пусть значения приближаемой функции f(x) заданы в N+1 узлах f(x0), ..., f(xN). Аппроксимирующую функцию будем выбирать из некоторого параметрического семейства F(x, c), где c = (c0, ..., cn)T — вектор параметров, N > n.

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

В этом случае задача аппроксимации ставится как задача поиска такого вектора параметров c = (c0, ..., cn)T, при котором значения аппроксимирующей функции как можно меньше отклонялись бы от значений аппроксимируемой функции  F(x, c) в совокупности всех узлов.

Графически задачу можно представить так


Запишем критерий среднеквадратичного приближения для метода наименьших квадратов:
J( c) = √ (Σi=0N[f(xi) — F(x, c) ]2) →min

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

Метод наименьших квадратов может быть применен для различных параметрических функций, но часто в инженерной практике в качестве аппроксимирующей функции используются многочлены по какому-либо линейно независимому базису {φk(x), k=0,...,n}:
F(x, c)  = Σk=0n[ckφk(x)].

В этом случае система линейных алгебраических уравнений для определения коэффициентов будет иметь вполне определенный вид:
a00c0 + a01c1 +… + a0ncn = b0
a10c0 + a11c1 +… + a1ncn = b1

an0c0 + an1c1 +… + anncn = bn

akj =  Σi=0N    [φk(xij(xi) ], bj =  Σi=0N[f(xij(xi) ]

Чтобы эта система имела единственное решение необходимо и достаточно, чтобы определитель матрицы А (определитель Грама) был отличен от нуля. Для того, чтобы система имела единственное решение необходимо и достаточно чтобы система базисных функций φk(x), k=0,...,n была линейно независимой на множестве узлов аппроксимации.

В этой статье рассматривается среднеквадратичное приближение многочленами по степенному базису  {φk(x) = xk, k=0,...,n}.

Пример


А теперь перейдем к примеру. Требуется  вывести эмпирическую формулу для приведенной табличной зависимости  f(х), используя метод наименьших квадратов.
x 0,75 1,50 2,25 3,00 3,75
y 2,50 1,20 1,12 2,25 4,28

Примем в качестве аппроксимирующей функцию
y = F(x) = c0 + c1x + c2x2, то есть, n=2, N=4

Система уравнений для определения коэффициентов:
a00c0 + a01c1 +… + a0ncn = b0
a10c0 + a11c1 +… + a1ncn = b1

an0c0 + an1c1 +… + anncn = bn

akj =  Σi=0Nk(xij(xi) ], bj =  Σi=0N[f(xij(xi) ]

Коэффициенты вычисляются по формулам:
a00 = N + 1 = 5, a01 = Σi=0Nxi = 11,25,  a02 = Σi=0Nxi2 = 30,94
a10 = Σi=0Nxi = 11,25,  a11 = Σi=0Nxi2 = 30,94,  a12 = Σi=0Nxi3 = 94,92
a20 = Σi=0Nxi2 = 30,94,  a21 = Σi=0Nxi3 = 94,92,  a22 = Σi=0Nxi4 = 303,76
b0 = Σi=0Nyi = 11,25,  b1 = Σi=0Nxiyi = 29,  b2 = Σi=0Nxi2yi = 90,21


Решаем систему уравнений и получаем такие значения коэффициентов:
c0 = 4,822, c1 = -3,882, c2 = 0,999

Таким образом
y = 4,8 — 3,9x + x2

График получившейся функции


Релизация на C#


А теперь перейдем к тому, как написать код, который бы строил такую матрицу. А тут, оказывается, все совсем просто:
private double[,] MakeSystem(double[,] xyTable, int basis)
{
  double[,] matrix = new double[basis, basis + 1];
  for (int i = 0; i < basis; i++)
  {
    for (int j = 0; j < basis; j++)
    {
      matrix[i, j] = 0;
    }
  }
  for (int i = 0; i < basis; i++)
  {
    for (int j = 0; j < basis; j++)
    {
      double sumA = 0, sumB = 0;
      for (int k = 0; k < xyTable.Length / 2; k++)
      {
        sumA += Math.Pow(xyTable[0, k], i) * Math.Pow(xyTable[0, k], j);
        sumB += xyTable[1, k] * Math.Pow(xyTable[0, k], i);
      }
      matrix[i, j] = sumA;
      matrix[i, basis] = sumB;
    }
  }
  return matrix;
}

На входе функция получает таблицу значений функций — матрицу, в первом столбце которой содержатся значения x, во втором, соответственно, y, а также значение степенного базиса.

Сначала выделяется память под матрицу, в которую будут записаны коэффициенты для решения системы линейных уравнений. Затем, собственно, составляем матрицу — в sumA записываются значения коэффициентов aij, в sumB — bi, все по формуле, указанной выше в теоретической части.

Для решения составленной системы линейных алгебраических уравнений в моей программе используется метод Гаусса. Архив с проектом можно скачать по ссылке.

Скриншот работы программы на примере, решенном выше:



Используемые источники:
Сулимова В.В. Методические указания по курсу «Вычислительный практикум» — Тула, ТулГУ, 2009 — 65 с.
Поделиться публикацией
Реклама помогает поддерживать и развивать наши сервисы

Подробнее
Реклама
Комментарии 11
  • +2
    График!!!
    • 0
      Спасибо за конструктивное замечание, добавил
    • 0
      2 вопроса:
      1 что такое dataGridTable?
      2 зачем заполнять масив нулями?
      • 0
        1. dataGridTable — таблица, в которую вводятся значения функции. Исправил «dataGridTable.RowCount» на «xyTable.Length / 2».
        2. В принципе да, абсолютно не обязательно, он итак полностью заполнится. Но с инициализацией недежнее.
      • +2
        Интересные здесь люди: то за «почти» велосипед минусуют, как только могут, то за явный — плюсуют.

        Между прочим, приближение функций методом наименьших квадратов элементарно (для дифференцируемых функций, конечно) реализуется библиотекой GSL. Автор, судя по всему, о ней не слышал.
        • –1
          Только вот GSL на С (если вы конечно про это: www.gnu.org/s/gsl/), а примеры тут на C#, и есть вероятность, что автору язык C не подходит…
          • +4
            Вот написал бы он на PHP — из минусов бы не вылез))

            На самом деле, вдимо не все IT-специалисты с хабра знакомы с численными методами. Квадратичная аппроксимация с точки зрения сложности реализации алгоритма близка сортировке пузырьком… (имхо)
            • –1
              Согласен, метод довольно простой. Но с другой стороны, это мой первый топик на Хабрахабре, и я решил начать с чего-то не очень сложного.
              • +1
                Вам не нужно оправдываться. Сам факт того, то вы взяли и написали статью характеризует вас только с лучшей стороны.

                Больше удивляет то, что если человек пишет статью про систему кэширования на PHP, тут же налетает толпа народа, яростно ставят минусы и пишут гадости в комментариях.

                А вообще, я за популяризацию математики и мат.методов в частности, так что если вы продолжитецикл статей, будет просто замечательно.
          • 0
            На самом деле, здесь совершенно не важно, какой у нас базис, от скольких переменных функции базиса, непрерывны они или нет… Достаточно того, что в точках X, в которых нам известно значение аппроксимируемой функции Y, мы знаем значения y_i=f_i(X) (f_i — базисные функции). Дальше с каждой точкой делаем как в статье — M[i,j]+=y_i*y_j для 0<=i,j<N (N — размерность базиса), M[i,N]+=y_i*Y. Решаем систему и получаем коэффициенты.
            К сожалению, мимнимизация суммы квадратов хорошо работает только для нормально распределенной ошибки. А в реальной жизни они не такие — всегда встречаются какие-нибудь выбросы (имеющие для данного метода разрушительные последствия), или наоборот, функция плотности оказывается многогорбой. Может быть, кто-нибудь знает хороший метод, как минимизировать, скажем, медиану значений abs(sum(a-i*f_i(X))-Y)?
            • 0
              Через матрицы предложенный алгоритм можно в три строчки поместить
              (иногда так понятнее, что происходит, да и в MathCAD проще решить/проверить алгоритм):

              1) Исходная СЛАУ (система алг. лин. уравненй) в матричном виде A*С=B — получена из таблицы опорных точек.

              2) Домножим слева на транспонированную A --> (AT*A)*С=(AT*B) получена СЛАУ, где AT*A=D — квадратная матрица размером N*N, а AT*B=F — вектор.

              3) Отсюда решение C=(AT*A)^(-1)*(AT*B), проще С=D^(-1)*F или Гаусс (если решается не в мат. пакете, а алгоритмически)

              Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.