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

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

Теория

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

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

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

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


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

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

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

В этом случае система линейных алгебраических уравнений для определения коэффициентов будет иметь вполне определенный вид:
a₀₀c₀ + a₀₁c₁ +… + a_0nc_n = b₀
a₁₀c₀ + a₁₁c₁ +… + a_1nc_n = b₁
…
a_n0c₀ + a_n1c₁ +… + a_nnc_n = b_n

a_kj =  Σ_i=0^N    [φ_k(x_i)φ_j(x_i) ], b_j =  Σ_i=0^N[f(x_i)φ_j(x_i) ]

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

В этой статье рассматривается среднеквадратичное приближение многочленами по степенному базису  {φ_k(x) = x^k, 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) = c₀ + c₁x + c₂x², то есть, n=2, N=4

Система уравнений для определения коэффициентов:
a₀₀c₀ + a₀₁c₁ +… + a_0nc_n = b₀
a₁₀c₀ + a₁₁c₁ +… + a_1nc_n = b₁
…
a_n0c₀ + a_n1c₁ +… + a_nnc_n = b_n

a_kj =  Σ_i=0^N[φ_k(x_i)φ_j(x_i) ], b_j =  Σ_i=0^N[f(x_i)φ_j(x_i) ]

Коэффициенты вычисляются по формулам:
a₀₀ = N + 1 = 5, a₀₁ = Σ_i=0^Nx_i = 11,25,  a₀₂ = Σ_i=0^Nx_i² = 30,94
a₁₀ = Σ_i=0^Nx_i = 11,25,  a₁₁ = Σ_i=0^Nx_i² = 30,94,  a₁₂ = Σ_i=0^Nx_i³ = 94,92
a₂₀ = Σ_i=0^Nx_i² = 30,94,  a₂₁ = Σ_i=0^Nx_i³ = 94,92,  a₂₂ = Σ_i=0^Nx_i⁴ = 303,76
b₀ = Σ_i=0^Ny_i = 11,25,  b₁ = Σ_i=0^Nx_iy_i = 29,  b₂ = Σ_i=0^Nx_i²y_i = 90,21

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

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

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

Релизация на 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 с.