Аппроксимация. Метод наименьших квадратов

Пусть данные некоторого эксперимента представлены в виде таблицы значений независимой переменной x и зависимой переменной y.

Требуется отыскать аналитическую зависимость f(x,a₀,a₁,…,a_m), являющуюся функцией одной независимой переменной x и параметров a₀,a₁,a₂,…,a_m, которая наилучшим образом описывала бы эти экспериментальные данные в смысле минимума квадратичного критерия рассогласования R(a₀,a₁,…,a_m):

Функцию f(x,a₀,a₁,…,a_m) определим как полином степени m вида: Надо найти такие значения параметров, при которых квадратичный критерий рассогласования имел бы минимальное значение  Вывод формулы для определения параметров в матричном виде рассмотрим на примере полинома второй степени (m=2).  Тогда критерий R будет являться функцией трёх переменных a₀, a₁, a₂:  Необходимые условия минимума критерия R имеют вид:

или

Полученную линейную относительно искомых параметров a0,a1,a2, систему уравнений запишем в матричном виде:  где

 Для удобства формирования матрицы коэффициентов  и столбца свободных членов  введем матрицу  элементы которой определяются через значения независимой переменной x_i, i=0,1,2,…,n  тогда

При аппроксимации полиномами высших порядков матрица  будет иметь вид:

В общем случае количество строк в матрице  равно количеству точек, а количество столбцов равно количеству параметров, где строка состоит из значений частных производных от функции f(x,a₀,a₁,…,a_m) по соответствующему параметру.

 

 

Блок-схема аппроксимации.

Begin End m,n,x,y a,dy,r a=inv(F’*F)*F’*y i=1 шаг 1 до n j=1 шаг 1 до m F(i,j)=x(i)^(j-1) yr=F*a; dy=y-yr r=dy’*dy

Пример. xi= -1 0 1 2;  yi= 2 1 2 4.

 

 

Интегрирование

Пусть на отрезке [a; b] определена непрерывная функция f(x). f=inline(‘<функция>'); Требуется определить значение определенного интеграла  которое числено равно площади S фигуры, ограниченной графиком функции f(x) и осью x, на заданном отрезке [a; b]. Для приближенного вычисления площади, разобьем отрезок [a; b] на n равных элементарных отрезков точками: x₀=a, x₁= a+h, x₂=x₁+h,…,x_i=x_i–1+h,…,x_n=b, – шаг разбиения. Значение функции f(x) в точках (рис.2.6.1) разбиения x_i обозначим через y_i.

x0=a x1 x3 x2 xn=b f(x) xi xn–1 xi+1 x · · ·    · · ·    · · ·   · · ·   y0 y1 y2 y3 yi–1 yi yi+1 yn–1 yn s0 s1 s2 sn-1 si xi–1 si-1

Рис.2.6.1. Интегрирование. Разбиение на равные отрезки.

f=inline(‘<функция>');

x=a:h:b;

plot(x,f(x),'k-')

Площадь S можно вычислить как сумму элементарных площадей определенных для соответствующих элементарных отрезков длиной h: S = s₀+s₁+s₂+…s_i+…..+s_n–1 Произвольную площадь s_i можно вычислить, как определенный интеграл на отрезке [x_i;x_i+1] от более простой функции φ_i(x), которой заменим реальную функцию f(x):  Вид функции φ_i(x) будет определять название метода.

Методы прямоугольников. Значение функции φi(x) на отрезке [x_i;x_i+1] принимается константой

Метод прямоугольников вперед. Для функции φ_i(x) = y_i значения элементарной s_i и общей S площади можно вычислить как:  , тогда

x=a:h:b-h;

S=h*sum(f(x));

Метод прямоугольников назад. Для функции φ_i(x) = y_i значения элементарной s_i и общей S площади можно вычислить как:  , тогда

x=a+h:h:b;

S=h*sum(f(x));

Метод прямоугольников в среднем. Вычислим  и значение функции  Тогда значения элементарной s_i и общей S площади можно вычислить как

x=a+h/2:h:b;

S=h*sum(f(x));

Метод трапеций. Функцию φ_i(x) будем определять как линейную на отрезке [x_i;x_i+1], т.е. ее график должен проходить через две смежные точки (x_i,y_i) и (x_i+1,y_i+1). Функцию φ_i(x) можно будет представить как интерполяционный многочлен Лагранжа, построенный по двум точкам (x_i,y_i) и (x_i+1,y_i+1):  тогда значения элементарной s_i площади можно вычислить как:  Введем переменную  Тогда x = x_i + h·t и dx = h·dt. Значениям x, равным x, x_i+1 соответствуют значения t, равные 0, 1. Значение (x-x_i) = x_i–x_i + h·t = h·t. Значение (x-x_i+1) = x_i – x_i+1+ h·t = h(t-1). Элементарную площадь s_i с использованием новой переменной определим как:

x=a:h:b-h; S=h*sum((f(x)+f(x+h))/2);

x=a:h:b; S=trapz(x,f(x));

x=a:h:b; S=h*trapz(f(x));

Метод Симпсона. Определим точку x_i+½ = x_i+½·h в середине элементарного отрезка [x_i;x_i+1] и значение функции в этой точке y_i+½ Функцию φ_i(x) будем определять как квадратичную на отрезке [x_i;x_i+1], т.е. её график должен проходить через три смежные точки (x_i,y_i),(x_i+½ , y_i+½) и (x_i+1,y_i+1). Функцию φ_i(x) можно будет представить как интерполяционный многочлен Лагранжа, построенный по трём точкам x_i, x_i+½ и x_i+1:

 Тогда значения элементарной s_i площади можно вычислить как:

 

 Введем переменную  тогда x = x_i + h·t и dx = h·dt. Значениям x, равным x_i, x_i+½, x_i+1 соответствуют значения t, равные 0,½,1 Значение (x-x_i) = x_i–x_i + h·t = h·t. Значение (x-x_i+½) = x_i – x_i+½ + h·t = h(t- ½) Значение (x-x_i+1) = x_i – x_i+1+ h·t = h(t-1) Элементарную площадь s_i с использование новой переменной определим как:

Тогда значения общей S площади можно вычислить как:

x=a+h:h:b-h; xs=a+h/2:h:b; S=h/6*(f(a)+f(b)+2*sum(f(x))+4*sum(f(xs)));

S=quad(f,a,b);