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

Пусть данные некоторого эксперимента представлены в виде таблицы значений независимой переменной 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) по соответствующему параметру.

 

Программа. Аппроксимация.

 

function DATA

global x y m n;

x=[-4:4];

y=[254;75;12;0;1;2;21;92;287];

n=9;

m=4;

end

 

function [ a,dy,r,yr ] = fun_Approximation(x,y,m,n)

for i=1:n

for j=1:m+1

   F(i,j)=x(i)^(j-1);

end %for j

end %for i

Ft=F';

N=Ft*F;

b=Ft*y;

Nobr=N^(-1);

a=Nobr*b;

yr=F*a;

dy=y-yr;

r=dy'*dy;

 end % function

 

function GLAV_Approximation

global x y m n a dy r yr;

DATA;

[ a,dy,r,yr ] = fun_Approximation(x,y,m,n);

REPORT;

end

 

function REPORT

global x y m n a dy r yr;

disp('Approximation');

disp(' ');

disp('     i    x        y      yr    dy'); 

disp(' ________________________________________________________')

i=0;

for i=1:length(x)

xx=x(i);

yy=y(i);

yrr=yr(i);

dyy=dy(i);

disp(sprintf('%10.1f\t%10.5f\t %10.5f\t%10.5f\t %10.5f',i,xx,yy,yrr,dyy));

end %for

disp(' ');

disp('   i    a '); 

disp('  _______________')

i=0;

for i=1:length(a)

aa=a(i);

disp(sprintf('%10.1f\t%10.5f',i,aa));

end %for

plot(x,y,'r*',x,yr,'k-');

grid on;

xlabel('x *, xx --');

ylabel('y *, yr --');

title('Approximation');

end

 

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

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);

 

Программа. Интегрирование. Сравнение погрешностей методов при одинаковом числе шагов

function DATA

global a b n;

a=-1;

b=4;

n=50;

end

 

function [ fx ] = f(x)

fx=x.^2-5;

end

 

function [ Integr ] = Pramougoln_Vpered(a,b,n)

disp('Pramougoln_Vpered');

h=(b-a)/n;

for i=1:n+1

if i==1

   x(i)=a;

else

   x(i)=x(i-1)+h;

end%if

y(i)=f(x(i));

end%for i

sum=0;

for i=1:n

sum=sum+y(i);

end%for i

Integr=h*sum;

end

 

function [ Integr ] = Pramougoln_Nazad(a,b,n)

disp('Pramougoln_Nazad');

h=(b-a)/n;

for i=1:n+1

if i==1

   x(i)=a;

else

   x(i)=x(i-1)+h;

end%if

y(i)=f(x(i));

end%for i

sum=0;

for i=2:n+1

sum=sum+y(i);

end%for i

Integr=h*sum;

end

 

function [ Integr ] = Pramougoln_Sredn(a,b,n)

disp('Pramougoln_Sredn');

h=(b-a)/n;

for i=1:n

if i==1

   x(i)=a+h/2;

else

   x(i)=x(i-1)+h;

end%if

y(i)=f(x(i));

end%for i

sum=0;

for i=1:n

sum=sum+y(i);

end%for i

Integr=h*sum;

end

 

function [ Integr ] = Trapez(a,b,n)

disp('Trapez');

h=(b-a)/n;

for i=1:n+1

if i==1

   x(i)=a;

else

   x(i)=x(i-1)+h;

end%if

y(i)=f(x(i));

end%for i

sum=0;

sum=(y(1)+y(n+1))/2;

for i=2:n

sum=sum+y(i);

end%for i

Integr=h*sum;

end

 

function [ Integr ] = Symps(a,b,n)

disp('Symps');

h=(b-a)/n;

for i=1:n+1

if i==1

   x(i)=a;

else

   x(i)=x(i-1)+h;

end%if

y(i)=f(x(i));

end%for i

for i=1:n

if i==1

   xsr(i)=a+h/2;

else

   xsr(i)=xsr(i-1)+h;

end%if

ysr(i)=f(xsr(i));

end%for i

sum=0;

for i=1:n

sum=sum+y(i)+4*ysr(i)+y(i+1);

end%for i

Integr=h*sum/6;

end

 

function GLAV_Integr

global a b n Integral;

DATA;

[ Integral(1) ] = Pramougoln_Vpered(a,b,n);

[ Integral(2) ] = Pramougoln_Nazad(a,b,n);

[ Integral(3) ] = Pramougoln_Sredn(a,b,n);

[ Integral(4) ] = Trapez(a,b,n);

[ Integral(5) ] = Symps(a,b,n);

REPORT;

end

 

function REPORT

global a b n Integral;

disp('Arguments f(x)=x^2-5');

disp(['left board a = ' num2str(a,'%10.5f') ]);

disp(['right board b = ' num2str(b,'%10.5f') ]);

disp(['Number of steps n = ' num2str(n,'%10.5f') ]);

disp(['Pramougoln_Vpered = ' num2str(Integral(1),'%10.5f') ]);

disp(['Pramougoln_Nazad = ' num2str(Integral(2),'%10.5f') ]);

disp(['Pramougoln_Sredn = ' num2str(Integral(3),'%10.5f') ]);

disp(['Trapez = ' num2str(Integral(4),'%10.5f') ]);

disp(['Symps = ' num2str(Integral(5),'%10.5f') ]);

end