Министерство образования Российской Федерации МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ ТЕХНИЧЕСКИЙ УНИВЕРСИТЕТ им. Н.Э. БАУМАНА






Лабораторная работа №2

По дисциплине “Численные методы”

Тема Метод наименьших квадратов







Выполнил: Студент группы МТ3-33

Пухов Р.А.


Проверил: Доцент кафедры ФН-1

Семакин А.Н.






Москва, 2015




Пусть известны значения yi в узлах xi , i = 0, 1, 2, ..., n, (a0, a1, ..., am, x) - функция, зависящая от параметров a0, a1, ..., am. Рассмотрим функцию S:

S=.


Выберем параметры a0, a1, ..., am так, чтобы минимизировать S , т.е. сумму квадратов невязок (Рис.1). Отсюда получаем систему уравнений:

= 0, i = 0, 1,…,m.

Эту систему уравнений (часто нелинейную) можно решить методом Ньютона. Рассмотрим подробнее случай, когда функция (x) является многочленом степени m:


Рис.1


Условие (1) приводит к следующей СЛАУ:



Преобразуем ее к виду:


Введем коэффициенты:

, .


Получим систему линейных алгебраическиx уравнений:



Эта система решается методом Гаусса.


Код программы:

clc

% Поиск значений Xi, Yi:

B=[0,0,0,0,0,0,0,0,0,0;

0,0,0,0,0,0,0,0,0,0];

for i=0:10

%Поиск значений X в интервале от a до b с шагом h

X(i+1)=(-1+0.2*i);

%Запись полученных значений в первую стору массива B

B(1,i+1)=X(i+1);

%Расчет значений Y(x) по полученным значениям X и запись их во вторую стоку массива B

B(2,i+1)=exp(X(i+1))-exp(2*X(i+1));

end;

B


%Поиск коэффициентов b(pq)

C=[0,0,0,0,0,0,0];

for j=0:6

for i=1:11

C(j+1)=C(j+1)+B(1,i)^j;

end

end

%Поиск коэффициентов с(p)

D=[0,0,0,0];

for i=1:11

D(1)=D(1)+B(2,i);

D(2)=D(2)+B(2,i)*B(1,i);

D(3)=D(3)+B(2,i)*B(1,i)^2;

D(4)=D(4)+B(2,i)*B(1,i)^3;

end

%Создание матрицы для расчёта кубической функции

Q=[0,0,0,0,0;

0,0,0,0,0;

0,0,0,0,0;

0,0,0,0,0];

Q(1,1)=C(1);

Q(1,2)=C(2);

Q(2,1)=C(2);

Q(1,3)=C(3);

Q(2,2)=C(3);

Q(3,1)=C(3);

Q(1,4)=C(4);

Q(2,3)=C(4);

Q(3,2)=C(4);

Q(4,1)=C(4);

Q(2,4)=C(5);

Q(3,3)=C(5);

Q(4,2)=C(5);

Q(3,4)=C(6);

Q(4,3)=C(6);

Q(4,4)=C(7);

Q(1,5)=D(1);

Q(2,5)=D(2);

Q(3,5)=D(3);

Q(4,5)=D(4);

%Создание матрицы для расчёта линейной функции

L=[0,0,0;

0,0,0];

L(1,1)=Q(1,1);

L(1,2)=Q(1,2);

L(2,1)=Q(2,1);

L(2,2)=Q(2,2);

L(1,3)=D(1);

L(2,3)=D(2);

%Создание матрицы для расчёта параболической функции

M=[0,0,0,0;

0,0,0,0;

0,0,0,0];

M(1,1)=Q(1,1);

M(1,2)=Q(1,2);

M(2,1)=Q(2,1);

M(2,2)=Q(2,2);

M(1,3)=Q(1,3);

M(2,3)=Q(2,3);

M(3,1)=Q(3,1);

M(3,2)=Q(3,2);

M(3,3)=Q(3,3);

M(1,4)=D(1);

M(2,4)=D(2);

M(3,4)=D(3);


%Поиск коффициентов для уравнений линейной, параболической и кубической функции

for i=1:3

%Условие для поиска линейной функции

if i==1

A=L;

N=2;

X=[0,0];

end

%Условие для поиска параболической функции

if i==2

A=M;

N=3;

X=[0,0,0];

end

%Условие для поиска кубической функции

if i==3

A=Q;

N=4;

X=[0,0,0,0];

end

%Решение матрицы методом Гаусса

for i=1:N

amax=A(i,i);

imax=i;

for j=(i+1):N

if abs(A(j,i))>abs(amax)

amax=A(j,i);

imax=j;

end;

end;

for j=1:(N+1)

b=A(i,j);

A(i,j)=A(imax,j);

A(imax,j)=b;

end;

b=A(i,i);

for j=1:(N+1)

A(i,j)=A(i,j)/b;

end;

for k=(i+1):N

b=A(k,i);

for j=i:(N+1)

A(k,j)=A(k,j)-A(i,j)*b;

end;

end;

end;

X(N)=A(N,N+1);

for i=1:(N-1)

i1=N-i;

sum=0;

for j=(i1+1):N

sum=sum+X(j)*A(i1,j);

end;

X(i1)=A(i1,N+1)-sum;

end;

z=X;

if i==1

X1=X;

end

if i==2

X2=X;

end

if i==3

X3=X;

end

end


%Вывод уравнений линейной, параболической и кубической функции

disp(sprintf('Линейная функция: f(x)= %g + %g*x',X1(1),X1(2)))

disp(sprintf('Параболическая функция: f(x)= %g + %g*x + %g*x^2',X2(1),X2(2), X2(3)))

disp(sprintf('Кубическая функция: f(x)= %g + %g*x + %g*x^2 + %g*x^3',X3(1),X3(2), X3(3), X3(4)))



%Подсчет и вывод ответа

O=[0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0;0,0,0,0,0,0,0,0,0];

for i=1:11

O(i,1)=B(1,i);

O(i,2)=B(2,i);

O(i,3)=X1(1)+X1(2)*B(1,i);

O(i,4)=X2(1)+X2(2)*B(1,i)+X2(3)*B(1,i)^2;

O(i,5)=X3(1)+X3(2)*B(1,i)+X3(3)*B(1,i)^2+X3(4)*B(1,i)^3;

O(i,6)=(X1(1)+X1(2)*B(1,i))-B(2,i);

O(i,7)=(X2(1)+X2(2)*B(1,i)+X2(3)*B(1,i)^2)-B(2,i);

O(i,8)=(X3(1)+X3(2)*B(1,i)+X3(3)*B(1,i)^2+X3(4)*B(1,i)^3)-B(2,i);

O(i,9)=O(i,6)+O(i,7)+O(i,8);

end

O

Ответ к лабораторной работе:



















Чтобы не видеть здесь видео-рекламу достаточно стать зарегистрированным пользователем.
Чтобы не видеть никакую рекламу на сайте, нужно стать VIP-пользователем.
Это можно сделать совершенно бесплатно. Читайте подробности тут.