Интерполяция подразумевает, что интерполянт точно проходит через заданные точки. Однако зачастую, особенно при обработке экспериментальных данных, когда сами узловые значения известны лишь приближённо, более разумным оказывается строить плавную кривую, неким оптимальным образом проходящую вблизи заданных точек, а не точно через них. Поиск такой функциональной зависимости называют "сглаживанием" экспериментальных данных.
Задачу о сглаживании экспериментальных данных можно решать, используя метод наименьших квадратов. Согласно методу наименьших квадратов указывается вид эмпирической формулы
где 
- числовые параметры.
Наилучшими значениями параметров считаются те, для которых сумма квадратов уклонений функции 
от экспериментальных точек (x i, y i) (i=1,2,...,m) является минимальной.
Если искомая сглаживающая функция – полином, то поиск его коэффициентов в MATLAB можно осуществить с помощью функции polyfit. Рассмотрим следующий пример.
time = zeros(0,2);
n = 100000; h = 5000; N = 2000000;
for m = n:h:N,
m,
a = rand(1,m);
tic;
sort(a);
time = [time; m, toc];
End
Save -ascii sorttime.dat time
Load sorttime.dat
time = sorttime;
x = time(:,1); y = time(:,2);
% x * log x
p = polyfit(x.* log(x), y, 1);
yy = polyval(p, x.* log(x));
err = sqrt(sum((yy - time(:, 2)).^2))
% This is the same as F.normr where F is such that
% [p, F] = polyfit(x.* log(x), y, 1);
Plot(x, y, '.', x, yy, 'r', 'LineWidth',2)
title('x log x');
Pause
% Further...
% linear function
[p, F] = polyfit(x, y, 1)
yy = polyval(p, x);
Plot(x, y, '.', x, yy, 'm', 'LineWidth',2)
title('ax+b')
Pause
% quadratic function
[p, F] = polyfit(x, y, 2)
yy = polyval(p, x);
Plot(x, y, '.', x, yy, 'b', 'LineWidth',2)
title('ax^2+bx+c')
Pause
% ax^b
[p, F] = polyfit(log(x), log(y), 1)
yy = exp(polyval(p, log(x)));
Plot(x, y, '.', x, yy, 'k', 'LineWidth',2)
title('ax^b')
Pause
% В следующих задачах нужно проводить отбор зависимости с
%минимальной степенью многочлена, для которой
% выполнено неравенство невязка < delta, где delta -
% среднеквадратичная относительная ошибка измерений.
%
% Зависимость модуля деформации материала у от нагрузки х.
% Данные получены с (относительной) точностью 1% -
% (delta=0.01).
x=[50 70 90 100 200]'; y=[51.546 61.838 73.230 77.339 115.607]';
t=[50:200]';
% Найдите формулы и графики этой зависимости, используя
% МНК для разных степеней полинома. Какую зависимость
% следует выбрать и почему?
% Какая из зависимостей не физична?
% Экспериментальная зависимость излучательной способности
% вольфрама у от длины волны излучения х. Данные получены
% с точностью 0.1% (delta=0.001).
х = [0.55 0.575 0.6 0.625 0.65 0.656 0.675 0.7 0.725 0.75 0.8 0.85 0.9]';
у = [0.45 0.445 0.441 0.438 0.434 0.433 0.43 0.426 0.422...
0.418 0.408 0.399 0.39]';
t=[0.55:0.005:0.9]';
% Найдите формулы и графики этой зависимости, используя
% МНК для разных степеней полинома. Какую зависимость
% следует выбрать и почему?
Вместо функции polyfit можно непосредственно составлять матрицу системы (матрица Вандермонда), которая будет переопределённой, и решать её с помощью оператора \ или, что одно и то же, рассчитав псевдообратную матрицу системы (функция pinv), как в примере ниже.
% Случай, когда полиномиальная зависимость для
% приближения не годится!
%
x=[0:12]'; y=10.^x+0.001*randn(size(x));
% Приближение полиномом степени 4:
A=[ones(size(x)) x x.^2 x.^3 x.^4];
% Степень полинома - макс, степень в матрице.
a=pinv(A)*y % найденные коэффициенты
yy=A*a; conl=cond (A)
nev=norm(y-yy)/norm(y) % невязка плохая!
%
% Сделаем по-другому: прологарифмируем данные у
Y=log10(y);
% Теперь применим МНК
aa=pinv(A)*Y % найденные коэффициенты
YY=A*aa; yyl=10.^YY;% Получили зависимость вила y=10^{ax^4+bxA3+... }
nev=norm(y-yyl)/norm(y) % Посмотрите невязку!
figure(1); subplot (2,1,1);plot (x,y, '.-', x, yy, 'or');
legend('Эксп. данные','y=ax^4+bx^3+...', 2);
subplot(2,1,2); plot(x,y, '.-', x,yyl,'or');
legend('Эксп. данные', 'y=10^{ax^4+bx^3+...}', 2);
%
% Обратите внимание на то, что первый метод дает значения yy < 0!
% Например,
yy(9)
% В зависимости вида у=10^х это абсурдно.
Литература.
1. Калиткин Н.Н., Численные методы. – М.: Наука, 1978.
2. Золотых Н.Ю. Matlab в научных исследованиях, 2004.
Задание.
1. Для функции Рунге R (x)=1/(1+25 x 2) на отрезке [-1,1] найти коэффициенты интерполяционного многочлена степени n по узлам, распределённым a) равномерно, б) по формуле Чебышева. В одних осях по 100 точкам построить графики исходной функции и найденных полиномов, отобразив также маркерами соответствующих цветов узлы интерполяции. На другом графике изобразить ошибки интерполяции для случаев а) и б). Вывести значения максимальных по модулю ошибок.
2. Проведя вычисления для диапазона n = 4:2:16, найти зависимость максимальной ошибки интерполяции от порядка интерполирующего полинома (количества узлов) в виде графика или таблицы для случаев а) и б);
3. Проанализировать результаты.
4. Пусть приближаемая функция – полином степени m. Как будут себя вести ошибки интерполяции полиномами степени n при n<m, n=m и n>m?
5. Выполнить упражнение с данными из следующей таблицы
| x | ||||||
| y1 | 4,8 | 5,76 | 6,912 | 8,294 | 9,95 | |
| y2 | 3,08 | 4,3 | 5,16 | 5,83 | ||
| y3 | 0,33 | 0,5 | 0,6 | 0,67 | 0,71 | |
| y4 | 1.5 | 1,75 | 1,83 | 1,87 | 1.9 | |
| y5 | 0,2 | 0,11 | 0,077 | 0,059 | ||
| y6 | 0.4 | 0,33 | 0,31 | 0,29 | ||
| y7 | 2,25 | 3.37 | 5,06 | 7,59 | 11.4 | |
| y8 | 2,69 | 3,1 | 3,39 | 3,61 |
1. Нанести экспериментальные точки (xi,yi) на координатную сетку (х,у).
2. Выбрать одну из шести предложенных формул преобразования к переменным (х,у) так, чтобы преобразованные экспериментальные данные (xi,yi) наименее уклонялись от прямой.
3. Методом наименьших квадратов найти наилучшие значения параметров k и b в уравнении прямой y=kx+b.
4. Найти явный вид эмпирической формулы y=Q(x,k,b) и построить график эмпирической функции.
