Компьютерное моделирование при обработке опытных данных

N=n+1.

Построение интерполяционного многочлена в явном виде

Для построения интерполяционного многочлена вида (11.2) необходимо определить его коэффициенты a0, a1, :, an, т.е. ai i=0,1,2,:,n. Количество неизвестных коэффициентов равно

n+1=N,

где

n-степень многочлена (11.2),

N-количество узловых точек табличной функции (11.1).

Для нахождения коэффициентов, используем свойство (11.3) интерполяционного многочлена (11.2). На основании этого свойства интерполяционный многочлен должен пройти через каждую узловую точку (xi, yi) таблицы (11.1), т.е.,

$$a_0x_i^n + a_1x_i^{n-1} + \ldots + a_{n-1}x_i + a_n = y_i, i=0,1,\ldots,n.$$

Подставляя в (11.4) каждую узловую точку таблицы (11.1) получаем систему линейных уравнений:

$$\left\{ \begin{array}{l} a_0x_0^n + a_1x_0^{n-1} + \ldots + a_{n-1}x_0 + a_n = y_0,\\ a_0x_1^n + a_1x_1^{n-1} + \ldots + a_{n-1}x_1 + a_n = y_1,\\ \ldots \ldots \ldots \ldots \ldots \ldots\\ a_0x_n^n + a_1x_n^{n-1} + \ldots + a_{n-1}x_n + a_n = y_n. \end{array} \right.$$

Неизвестными системы (11.5) являются a0, a1, a2, :, an т.е. коэффициенты многочлена (11.2). Коэффициенты при неизвестных системы (11.5) $$x_i^n, x_i^{n-1}, \ldots, x_i^0, i=0,1, \ldots,n, i=0,1,:,n$$ легко могут быть определены на основании данных таблицы (11.1).

Интерполяция по Лагранжу

Интерполяционный многочлен может быть построен при помощи специальных интерполяционных формул Лагранжа, Ньютона, Стерлинга, Бесселя и др.

Интерполяционный многочлен по формуле Лагранжа имеет вид:

$$L_n(x)=\frac{(x-x_1)(x-x_2)(x-x_3) \ldots \ldots (x-x_n)}{(x_0-x_1)(x_0-x_2)(x_0-x_3) \ldots (x_0-x_n)} \cdot y_0 +\\ + \frac{(x-x_0)(x-x_2)(x-x_3) \ldots \ldots (x-x_n)}{(x_1-x_0)(x_1-x_2)(x_1-x_3) \ldots (x_1-x_n)} \cdot y_1 +\\ + \frac{(x-x_0)(x-x_1)(x-x_3) \ldots \ldots (x-x_n)}{(x_2-x_0)(x_2-x_1)(x_2-x_3) \ldots (x_2-x_n)} \cdot y_2 + \ldots\\ + \frac{(x-x_0)(x-x_1)(x-x_1) \ldots \ldots (x-x_{n-1})}{(x_n-x_0)(x_n-x_1)(x_n-x_1) \ldots (x_n-x_{n-1})} \cdot y_n.$$

Докажем, что многочлен Лагранжа является интерполяционным многочленом, проходящим через все узловые точки, т.е. в узлах интерполирования xi выполняется условие Ln(xi) = yi. Для этого будем последовательно подставлять значения координат узловых точек таблицы (11.1) в многочлен (11.6). В результате получим:

если x=x0, то Ln(x0) = y0,

если x=x1, то Ln(x1) = y1,

:::::

если x=xn, то Ln(xn) = yn.

Это достигнуто за счет того, что в числителе каждой дроби при соответствующем значении уj, j=0,1,2,:,n отсутствует сомножитель (x-xi), в котором i=j, а знаменатель каждой дроби получен заменой переменной х на соответствующее значение хj.

Таким образом, интерполяционный многочлен Лагранжа приближает заданную табличную функцию, т.е. Ln(xi) = yi и мы можем использовать его в качестве вспомогательной функции для решения задач интерполирования, т.е. $$L_n(x_k) \approx y_k$$.

Чем больше узлов интерполирования на отрезке [x0,xn], тем точнее интерполяционный многочлен приближает заданную табличную функцию (11.1), т.е. тем точнее равенство:

$$f(x_k) \approx L_n(x_k).$$

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

При решении задачи экстраполирования функции с помощью интерполяционного многочлена вычисление значения функции за пределами отрезка [x0,xn] обычно производят не далее, чем на один шаг h, равный наименьшей величине$$\left|x_{i+1} - x_i\right|,$$ так как за пределами отрезка [x0,xn] погрешности, как правило, увеличиваются.

Программирование формулы Лагранжа

Свернем формулу Лагранжа (11.6). В результате получим

$$Ln(x) = \sum \limits_{j=0}^{n} B_j \cdot y_j,$$

где

$$B_j = \prod \limits_{i=0}^{n} \frac{x-x_i}{x_j-x_i},$$

но при этом обязательно выполнение условия $$i \neq j$$.

При построении алгоритма используют конструкцию из двух включенных циклов:

Внешним циклом накапливаем сумму $$L = \sum \limits_{j=0}^{n} B_j \cdot y_j$$.

Внутренним циклом накапливаем произведение $$B_j = \prod \limits_{i=0}^{n} \frac{x-x_i}{x_j-x_i}, i \neq j$$.

Алгоритм (рис.11.3) не предусматривает получение интерполяционного многочлена в явном виде, а сразу решает задачу интерполирования функции в заданной точке, x=D.

Обозначения в алгоритме:

n - степень интерполяционного многочлена Лагранжа (11.6), равная количеству узловых точек N минус один, т.е. n=N-1.

D - значение аргумента в точке, для которой решается задача интерполирования табличной функции (11.1).

L - значение многочлена (11.6).

(рис 11.3) Схема алгоритма интерполяции по Лагранжу

Интерполяция по Ньютону

Дана табличная функция:

i	xi	yi
0	x0	y0
1	x1	y1
2	x2	y2
...	...	...
n	xn	yn

или

$$y_i = f(x_i), i=\overline{0,n}.$$

Точки с координатами (xi, yi) называются узловыми точками или узлами.

Количество узлов в табличной функции равно

N=n+1.

Необходимо найти значение этой функции в промежуточной точке, например, x=D, причем $$D \in[x_0,x_n]$$.

Для решения задачи строим интерполяционный многочлен.

Интерполяционный многочлен по формуле Ньютона имеет вид:

$$L_n(x) = f(x_0) + (x - x_0) \cdot f(x_0,x_1) +\\ + (x - x_0) \cdot(x - x_1) \cdot f(x_0,x_1,x_2) +\\ + (x - x_0) \cdot(x - x_1) \cdot(x - x_2) \cdot f(x_0,x_1,x_2,x_3) + \ldots +\\ + (x - x_0) \cdot(x - x_1) \cdot \ldots \cdot (x - x_{n-1}) \cdot f(x_0,x_1, \ldots,x_n),$$

где

n - степень многочлена,

$$f(x_0), f(x_0,x_1), f(x_0,x_1,x_2), f(x_0,x_1,\ldots, x_n)$$ - разделенные разности 0-го, 1-го, 2-го,:., n-го порядка, соответственно.

Разделенные разности

Значения f(x0), f(x1), : , f(xn), т.е. значения табличной функции в узлах, называются разделенными разностями нулевого порядка (k=0).

Отношение $$f(x_0,x_1) = \frac{f(x_1)-f(x_0)}{x_1 - x_0}$$ называется разделенной разностью первого порядка (k=1) на участке [x0, x1] и равно разности разделенных разностей нулевого порядка на концах участка [x0, x1], разделенной на длину этого участка.

Для произвольного участка [xi, xi+1] разделенная разность первого порядка (k=1) равна

$$f(x_i,x_{i+1}) = \frac{f(x_{i+1})-f(x_i)}{x_{i+1} - x_i}.$$

Отношение $$f(x_0,x_1,x_2) = \frac{f(x_1,x_2)-f(x_0,x_1)}{x_2 - x_0}$$ называется разделенной разностью второго порядка (k=2) на участке [x0, x2] и равно разности разделенных разностей первого порядка, разделенной на длину участка [x0, x2].

Для произвольного участка [xi, xi+2] разделенная разность второго порядка (k=2) равна

$$f(x_i,x_{i+1},x_{i+2}) = \frac{f(x_{i+1},x_{i+2}) - f(x_i,x_{i+1})}{x_{i+2} - x_i}$$

Таким образом, разделенная разность k -го порядка на участке [xi, xi+k] может быть определена через разделенные разности (k-1) -го порядка по рекуррентной формуле:

$$f(x_i,x_{i+1},x_{i+2}, \ldots,x_{i+k}) = \frac{ f(x_{i+1},x_{i+2}, \ldots,x_{i+k}) - f(x_i,x_{i+1}, \ldots,x_{i+k-1})}{x_{i+k} - x_i}$$

где

$$k=\overline{1,n},$$ $$i=\overline{0,n-k},$$

n - степень многочлена.

Максимальное значение k равно n. Тогда i =0 и разделенная разность n -го порядка на участке [x0,xn] равна $$f(x_0,x_i, \ldots,x_n) = \frac{ f(x_1,x_2, \ldots,x_n) - f(x_0,x_1, \ldots,x_n-1)}{x_n - x_0}$$, т.е. равна разности разделенных разностей (n-1) -го порядка, разделенной на длину участка [x0,xn].

Разделенные разности $$f(x_0,x_1), f(x_0,x_1,x_2), \ldots, f(x_0,x_1,\ldots, x_n)$$ являются вполне определенными числами, поэтому выражение (11.7) действительно является алгебраическим многочленом n -й степени. При этом в многочлене (11.7) все разделенные разности определены для участков [x0, x0+k], $$k=\overline{1,n}$$.

Лемма: алгебраический многочлен (11.7), построенный по формулам Ньютона, действительно является интерполяционным многочленом, т.е. значение многочлена в узловых точках равно значению табличной функции

$$L_n(x_i) = f(x-i) = y_i; i=0,1, \ldots n.$$

Докажем это. Пусть х=х0, тогда многочлен (11.7) равен

$$L_n(x_0) = f(x_0) = y_0.$$

Пусть х=х1, тогда многочлен (11.7) равен

$$L_n(x_1) = f(x_0) + (x_1 - x_0) \cdot f(x_0,x_1) = f(x_0) + (x_1 - x_0) \cdot \frac{f(x_1) - f(x_0)}{x_1 - x_0} = f(x_1) = y_1$$

Пусть х=х2, тогда многочлен (11.7) равен

$$L_n(x_2) = f(x_0) + (x_2 - x_0) \cdot f(x_0,x_1) + (x_2 - x_0) \cdot (x_2 - x_1) \cdot f(x_0,x_1,x_2) =\\= f(x_0) + (x_2 - x_0) \cdot \frac{f(x_1) - f(x_0)}{x_1 - x_0}+ (x_2 - x_0) \cdot (x_2 - x_1) \cdot \frac{f(x_1,x_2) - f(x_0,x_1)}{x_2 - x_0}= f(x_1) = y_1$$

Заметим, что решение задачи интерполяции по Ньютону имеет некоторые преимущества по сравнению с решением задачи интерполяции по Лагранжу. Каждое слагаемое интерполяционного многочлена Лагранжа зависит от всех значений табличной функции yi, i=0,1,:n. Поэтому при изменении количества узловых точек N и степени многочлена n (n=N-1) интерполяционный многочлен Лагранжа требуется строить заново. В многочлене Ньютона при изменении количества узловых точек N и степени многочлена n требуется только добавить или отбросить соответствующее число стандартных слагаемых в формуле Ньютона (11.7). Это удобно на практике и ускоряет процесс вычислений.

Программирование формулы Ньютона

Для построения многочлена Ньютона по формуле (11.7) организуем циклический вычислительный процесс по $$k=\overline{1,n}$$. При этом на каждом шаге поиска находим разделенные разности k -го порядка. Будем помещать разделенные разности на каждом шаге в массив Y.

Тогда рекуррентная формула (11.8) будет иметь вид:

$$y_i = \frac{y_{i+1} - y_i}{x_{i+k} - x_i}\\ k=\overline{1,n};\\ i=\overline{0,n-k}.$$

В формуле Ньютона (11.7) используются разделенные разности k -го порядка, подсчитанные только для участков [x0, x0+k], т.е. разделенные разности k -го порядка для i=0. Обозначим эти разделенные разности k-го порядка как у0. А разделенные разности, подсчитанные для I > 0, используются для расчетов разделенных разностей более высоких порядков.

Используя (11.9), свернем формулу (11.7). В результате получим

$$L_n(x) = y_o + \sum \limits_{k=1}^{n}P \cdot y_0^*$$

где

у0 - значение табличной функции (11.1) для x=x0.

$$y_0^*$$ - разделенная разность k-го порядка для участка [x0, x0+k].

$$P = (x - x_0)(x - x_1)\ldots (x - x_{k-1}) = \prod \limits_{j=0}^{k-1} (x - x_j).$$

Для вычисления Р удобно использовать рекуррентную формулу P = P(x - xk-1) внутри цикла по k.

Схема алгоритма интерполяции по Ньютону представлена на рис.11.4.

(рис 11.4) Схема алгоритма интерполяции по Ньютону

Пример интерполяции по Ньютону

Дана табличная функция:

i	xi	yi
0	2	0,693147
1	3	1,098613
2	4	1,386295
3	5	1,609438

Вычислить разделенные разности 1-го, 2-го, 3-го порядков (n=3) и занести их в диагональную таблицу.

Разделенные разности первого порядка:

$$f(x_0,x_1) = \frac{f(x_1) - f(x_0)}{x_1 - x_0}= \frac{1,098613 - 0,693147}{3 - 2} = 0,405466.\\ f(x_1,x_2) = \frac{f(x_2) - f(x_1)}{x_2 - x_1}= \frac{1,386295 - 1,098613}{4 - 3} = 0,287682.\\ f(x_2,x_3) = \frac{f(x_3) - f(x_2)}{x_3 - x_2}= \frac{1,609438 - 1,386295}{5 - 4} = 0,223143.$$

Разделенные разности второго порядка:

$$f(x_0,x_1,x_2) = \frac{f(x_1,x_2) - f(x_0,x_1)}{x_2 - x_0}= \frac {0,287682 - 0,405466}{4 - 2} = -0,058892.\\ f(x_1,x_2,x_3) = \frac{f(x_2,x_3) - f(x_0,x_2)}{x_3 - x_1}= \frac {0,223143 - 0,287682}{5 - 3} = - 0,0322695.$$

Разделенная разность третьего порядка:

$$f(x_0,x_1,x_2,x_3) = \frac{f(x_1,x_3) - f(x_0,x_2)}{x_3 - x_0}= \frac {- 0,0322695 - (- 0,058892)}{5 - 2} = 0,00887416$$

Диагональная таблица разделенных разностей

i	xi	Разделенная разность
		0-го пор.	1-го пор.	2-го пор.	3-го пор.
0	2	0,693147
			0,405466
1	3	1,098613		-0,058892
			0,287682		0,00887416
2	4	1,386295		-0,0322695
			0,223143
3	5	1,60943

Интерполяционный многочлен Ньютона для заданной табличной функции имеет вид:

$$L_3(x) = f(x_0) + (x - x_0) \cdot f(x_0,x_1) + (x-x_0)(x - x_1) \cdot f(x_0,x_1,x_2) +\\+ (x - x_0)(x - x_1)(x - x_2) \cdot f(x_0,x_1,x_2,x_3) = 0,693147 + (x - 2) \cdot 0,405466+\\ + (x-2)(x-3) \cdot (-0,058892) + (x-2)(x-3)(x-4) \cdot 0,0887416.$$

Далее полученный интерполяционный многочлен Ньютона можно привести к нормальному виду$$L_3(x) = a_0x^3 + a_1x^2 + a_2x + a_3$$ и использовать его для решения задач интерполирования или прогноза.

Сплайн-интерполяция

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

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

Например, для некоторых функций (рис.11.5) необходимо задать все кубические функции q1(x), q2(x), :qn(x).

В наиболее общем случае эти многочлены имеют вид:

$$q_i(x) =k_{1j} + k_{2i} x + k_{3i} x^2 + k_{4i} x^3, i=\overline{1,n},$$

где kij - коэффициенты, определяемые описанными ранее условиями, количество которых равно 4n. Для определения коэффициентов kij необходимо построить и решить систему порядка 4n.

(рис 11.5)

Первые 2n условий требуют, чтобы сплайны соприкасались в заданных точках:

$$q_i(x_i) = y_i, i=\overline{1,n};\\ q_{i+1}(x_i) = y_i, i=\overline{0,n-1}.$$

Следующие (2п-2) условий требуют, чтобы в местах соприкосновения сплайнов были равны первые и вторые производные:

$$q_{i+1}'(x_i) = q'_i(x_i), i=\overline{1,n-1};\\ q''_{i+1}(x_i) = q_i''(x_i), i=\overline{0,n-1}.$$

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

$$q_1''(x_0) = 0, ..., q_n''(x_n) = 0.$$

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

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

Аппроксимация опытных данных

В результате проведения натурного эксперимента получена табличная функция:

i	X	Y
0	xo	yo
1	x1	y1
2	x2	y2
3	x3	y3
:	:	:
n	xn	yn

где

N-количество узловых точек в таблице,

n=N-1.

Задача аппроксимации заключается в отыскании аналитической зависимости y=f(x) полученной табличной функции.

В настоящее время существует 2 способа аппроксимации опытных данных:

Первый способ. Этот способ требует, чтобы аппроксимирующая кривая F(x), аналитический вид которой необходимо найти, проходила через все узловые точки таблицы. Эту задачу можно решить с помощью построения интерполяционного многочлена степени n:

$$P_n(x) = a_0 x^n + a_1 x^{n-1} + a_2 x^{n-2} + \ldots + a_{n-1} x^1 + a_n.$$

Однако этот способ аппроксимации опытных данных имеет недостатки:

Точность аппроксимации гарантируется в небольшом интервале [x0, xn] при количестве узловых точек не более 7-8.

Значения табличной функции в узловых точках должны быть заданы с большой точностью.

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

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

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

Сглаживание опытных данных методом наименьших квадратов

В этом методе при сглаживании опытных данных аппроксимирующей кривую F(x) стремятся провести так, чтобы ее отклонения $$\varepsilon_i$$ от табличных данных (уклонения) по всем узловым точкам были минимальными (рис 11.6), т.е.

$$\varepsilon_i=\left|F(x_i) - y_i \right| \to min.$$ (рис 11.6)

Избавимся от знака уклонения. Тогда условие (11.6) будет иметь вид:

$$\varepsilon_i^2=\left|(F(x_i) - y_i)^2 \right| \to min.$$

Суть метода наименьших квадратов заключается в следующем: для табличных данных, полученных в результате эксперимента, отыскать аналитическую зависимость F(x), сумма квадратов уклонений которой от табличных данных по всем узловым точкам была бы минимальной, т.е.

$$\sum \limits_{i=0}^{n}\varepsilon_i^2 = \sum \limits_{i=0}^{n} (F(x_i) - y_i)^2 \to min$$

Для определенности задачи искомую функцию F(x) будем выбирать из класса алгебраических многочленов степени m:

$$P_m(x) = a_0 x^m + a_1 x^{m-1} + a_2 x^{m-2} + \ldots + a_{m-1} x^1 + a_m.$$

Назовем многочлен (11.9) аппроксимирующим многочленом. Аппроксимирующий многочлен не проходит через все узловые точки таблицы. Поэтому его степень m не зависит от числа узловых точек. При этом всегда m < n. Степень m может меняться в пределах $$1 \le m \le N-2$$.

Если m=1, то мы аппроксимируем табличную функцию прямой линией. Такая задача называется линейной регрессией.

Если m=2, то мы аппроксимируем табличную функцию квадратичной параболой. Такая задача называется квадратичной аппроксимацией.

Если m=3, то мы аппроксимируем табличную функцию кубической параболой. Такая задача называется кубической аппроксимацией.

Уточним метод наименьших квадратов: для табличной функции, полученной в результате эксперимента, построить аппроксимирующий многочлен (11.9) степени m, для которого сумма квадратов уклонений по всем узловым точкам минимальна, т.е.

$$S = \sum \limits_{0}^{n} (P_m(x_i) - y_i) \to min.$$

Изменим вид многочлена Pm. Поставим на последнее место слагаемые, содержащие xm. На предпоследнее - слагаемые, содержащие xm-1 и т.д. В результате получим:

$$P_m(x) = a_0 x^0+ a_1 x^1 + a_2 x^2 + \ldots + a_m x^m.$$

или

$$P_m(x) = \sum \limits_{j=0}^{m}a_j x^j$$

При этом изменим индексы коэффициентов многочлена. Тогда условие (11.8) будет иметь вид:

$$S = \sum \limits_{i=0}^{n} (a_0 x_i^0 + a_1 x_i^1 + a_2 x_i^2 + \ldots + a_m x_i^m - y_i)^2 \to min,$$

где

xi и yi - координаты узловых точек таблицы,

aj, $$j=\overline{0,m}$$ -неизвестные коэффициенты многочлена (11.11).

Необходимым условием существования минимума функции S является равенство нулю ее частных производных по каждой aj.

$$\left\{ \begin{array}{l} \frac{\delta s}{\delta a_0} \to 2\sum \limits_{i=0}^{n}((a_0 x_i^0 + a_1 x_i^1 + \ldots + a_m x_i^m - y_i)x_i^0) = 0\\ \frac{\delta s}{\delta a_1} \to 2\sum \limits_{i=0}^{n}((a_0 x_i^0 + a_1 x_i^1 + \ldots + a_m x_i^m - y_i)x_i^1) = 0\\ \frac{\delta s}{\delta a_2} \to 2\sum \limits_{i=0}^{n}((a_0 x_i^0 + a_1 x_i^1 + \ldots + a_m x_i^m - y_i)x_i^2) = 0\\ \ldots \ldots \ldots \ldots \ldots \ldots \ldots \ldots\\ \frac{\delta s}{\delta a_m} \to 2\sum \limits_{i=0}^{n}((a_0 x_i^0 + a_1 x_i^1 + \ldots + a_m x_i^m - y_i)x_i^m) = 0 \end{array} \right.$$

В результате получили систему линейных уравнений. Раскрывая скобки и перенося свободные члены в правой части уравнений, получим в нормальной форме систему линейных уравнений:

$$\left\{ \begin{array}{l} c_0 a_0 + c_1 a_1 + c_2 a_2 + \ldots + c_m a_m =d_0,\\ c_1 a_0 + c_2 a_1 + c_3 a_2 + \ldots + c_{m+1} a_m =d_1,\\ c_2 a_0 + c_3 a_1 + c_4 a_2 + \ldots + c_{m+2} a_m =d_0,\\ \ldots \ldots \ldots \ldots \ldots \ldots\\ c_m a_0 + c_{m+1} a_1 + c_{m+2} a_2 + \ldots + c_{2m} a_m =d_m, \end{array} \right.$$

где

aj - неизвестные системы линейных уравнений (11.12),

$$c_k = \sum \limits_{i=0}^{n} x_i^k, k=\overline{0,2m}$$ - коэффициенты системы линейных уравнений (11.12),

$$d_j = \sum \limits_{i=0}^{n} y_i x_i^k, j=\overline{0,m}$$ - свободные члены системы линейных уравнений (11.12),

Порядок системы равен m+1.

При ручном счете коэффициенты ck и свободные члены dj удобно определять, пользуясь таблицей 11.2:

i	xi0	xi1	xi2	...	xi2m	xi0 yi	xi1 yi	...	xim
0	1
1	1
2	1
...	...
N	1
$$\sum_{i=o}^n$$	c0	c1	c2	...	c2m	d0	d1	...	dm

Программирование метода наименьших квадратов (МНК)

Изменим индексацию в системе (11.12). В результате получим:

$$\left\{ \begin{array}{l} c_{11}a_1 + c_{12}a_2 + c_{13}a_3 + \ldots + c_{1(m+1)} a_{m+1} = d_1,\\ c_{21}a_1 + c_{22}a_2 + c_{23}a_3 + \ldots + c_{2(m+1)} a_{m+1} = d_2,\\ c_{31}a_1 + c_{32}a_2 + c_{33}a_3 + \ldots + c_{3(m+1)} a_{m+1} = d_3,\\ \ldots \ldots \ldots \ldots \ldots \ldots \ldots \ldots\\ c_{(m+1)},1a_1 + c_{(m+1)},2a_2 + \ldots + c_{(m+1),(m+1)} a_{m+1} = d_{m+1} \end{array} \right.$$

где

$$a_j, j=\overline{1,(m+1)}$$ - неизвестные системы линейных уравнений (11.13),

$$c_{k,j} = \sum \limits_{i=1}^{N} x_i^{k+j-2}, k = \overline{1,(m+1)}, j = \overline{1,(m+1)}$$ - коэффициенты системы линейных уравнений (11.13),

$$d_k = \sum \limits_{i=1}^{N} y_i x_i^{j-1}, j = \overline{1,(m+1)}$$ - свободные члены системы линейных уравнений (11.13),

(xi, yi) - координаты узловых точек табличной функции, $$i=\overline{1,N}$$,

N - количество узловых точек,

m - степень аппроксимирующего многочлена вида:

$$P_m(x) = a_1 x^0 + a_2 x^1 + a_3 x^2 + \ldots + a_{m+1} x^m.$$

Алгоритм задачи:

Строим систему линейных уравнений (11.13). Определяем коэффициенты ck,j и свободные члены dk. Т.к. система (11.13) симметрична относительно главной диагонали, то достаточно определить только наддиагональные элементы системы.

Решаем систему (11.13) методом Гаусса. Находим коэффициенты aj многочлена (11.14).

Строим аппроксимирующий многочлен (11.14) и определяем его значение в каждой узловой точке Pi = Pm(xi).

Находим уклонение каждой узловой точки $$\varepsilon_i = P_i - y_i$$.

Находим сумму квадратов уклонений по всем узловым точкам $$S = \sum \limits_{i=1}^{N} \varepsilon_i^2$$.

Находим остаточную дисперсию $$D=\frac{S}{N-(m+1)}$$.

Для построения аппроксимирующего многочлена (11.11) и вычисления его значения в каждой узловой точке используем рациональную форму многочлена:

$$P_m(x) = a_1 + x \cdot (a_2 + x \cdot (a_3 + \ldots + x \cdot (a_m + x \cdot a_{(m+1)} \ldots)$$

Тогда для вычисления значения многочлена (11.15) удобно пользоваться схемой Горнера. Рекуррентная формула по схеме Горнера имеет вид:

$$P = a_{m+1},\\ P = a_j + x_i P, j=\overline{m,1,-1}, i=\overline{1,N}.$$

Укрупненная схема алгоритма МНК представлена на рис.11.7. Схемы алгоритмов основных блоков представлены на рисунках 11.8-11.10.

(рис 11.7) Укрупненная схема алгоритма аппроксимации методом наименьших квадратов

Обозначения в блоке 2:

m - степень аппроксимирующего многочлена,

N - количество узловых точек таблицы (11.2),

X, Y - массивы значений x и y таблицы (11.2).

(рис 11.8) Схема алгоритма блока 3. Определение коэффициентов системы (11.13) (рис 11.9) Схема алгоритма блока 4. Определение свободных членов системы (11.13) (рис 11.10) Схема алгоритма блока 6. Схема Горнера