Аппроксимация МНК

Оглавление

Введение. 3

Теоретическая часть. 5

Постановка задачи. 5

Построение эмпирических формул методом наименьших квадратов. 6

Метод Гаусса решения систем линейных алгебраических уравнений. 9

Блок-схема алгоритма программы на языке Pascal 11

Практическая часть. 15

Текст программы на языке Pascal, перечень использованных в программе идентификаторов и полученные результаты.. 15

Описание решения задачи в среде MathCAD.. 19

Результаты вычислений и их анализ. 20

Заключение. 22

Список литературы.. 23

Введение

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

Как известно, между величинами может существовать точная (функциональная) связь, когда одному значению аргумента соответствует одно определенное значение.

При выборе аппроксимации следует исходить из конкретной задачи исследования. Обычно, чем более простое уравнение используется для аппроксимации, тем более приблизительно получаемое описание зависимости. Поэтому важно считывать, насколько существенны и чем обусловлены отклонения конкретных значений от получаемого тренда. При описании зависимости эмпирически определенных значений можно добиться и гораздо большей точности, используя какое-либо более сложное, многопараметрическое уравнение. Однако нет никакого смысла стремиться с максимальной точностью передать случайные отклонения величин в конкретных рядах эмпирических данных. Выбирая метод аппроксимации, исследователь всегда идет на компромисс: решает, в какой степени в данном случае целесообразно и уместно «пожертвовать» деталями и, соответственно, насколько обобщенно следует выразить зависимость сопоставляемых переменных. Наряду с выявлением закономерностей замаскированных случайными отклонениями эмпирических данных от общей закономерности, аппроксимация позволяет также решать много других важных задач: формализовать найденную зависимость; найти неизвестные значения зависимой переменной путем интерполяции или, если это допустимо, экстраполяции.

Целью данной курсовой работы является изучение теоретических основ аппроксимации табулированной функции методом наименьших квадратов, и, применяя теоретические знания, нахождение аппроксимирующих полиномов. Нахождение аппроксимирующих полиномов в рамках данной курсовой работы следует путем написания программы на языке Pascal, реализующую разработанный алгоритм нахождения коэффициентов аппроксимирующего полинома, а также решить эту же задачу средствами MathCad.

В данной курсовой работе программа на языке Pascal разработана в оболочке PascalABC версия 1.0 beta. Решение задачи в среде MathCad производили в Mathcad версия 14.0.0.163.

Теоретическая часть

Постановка задачи

В данной курсовой работе необходимо выполнить следующее:

1. Разработать алгоритм нахождения коэффициентов трёх аппроксимирующих полиномов (многочленов) вида для табулированной функции y=f(x):

x	0,15	-0,01	0,46	1,71	3,94	7,4	12,41	19,47	29,21	42,49
f(x)	7,26	7,95	8,77	9,72	10,78	12,45	14,74	18,18	23,36	30,99
x	60,43	84,42	116,18	157,8	211,78	281,04	369	479,59	617,3	787,21
f(x)	41,97	57,38	78,45	106,65	143,64	191,28	251,71	327,27	420,58	534,52
x	995,04	1247,17	1550,71	1913,5	2344,18	2852,21	3447,91	4142,52	4948,19	5878,09
f(x)	672,24	837,2	1033,14	1264,13	1534,54	1849,11	2212,91	2631,36	3110,26	3665,8
x	995,04	1247,17	1550,71	1913,5	2344,18	2852,21	3447,91	4142,52	4948,19	5878,09
f(x)	672,24	837,2	1033,14	1264,13	1534,54	1849,11	2212,91	2631,36	3110,26	3665,8
x	6946,36	8768,24	9560,03	11139,19	12924,34	14935,29	17193,14	19720,24	22540,29	25678,32
f(x)	4274,55	4973,5	5760,63	6641,98	7627,6	8725,62	9945,22	11296,05	12788,24	14432,44

для степени полиномов n=2, 4, 5.

2. Построить блок-схему алгоритма.

3. Создать программу на языке Pascal, реализующую разработанный алгоритм.

4. Рассчитать среднеквадратичные отклонения для каждого из трех случаев по формуле:

5. Построить графики 3-х полученных приближающих функций в одной системе координат. На графике должны содержаться и исходные точки (х
i
,
yi
)
.

6. Решить задачу средствами MathCAD.

Результаты решения задачи с помощью созданной программы на языке Pascal и в среде MathCAD нужно представить в виде построенных с помощью найденных коэффициентов трёх полиномов; таблицы, содержащей полученные с помощью найденных полиномов значения функции в точках хi и среднеквадратичных отклонений.

Построение эмпирических формул методом наименьших квадратов

Очень часто, особенно при анализе эмпирических данных возникает необходимость найти в явном виде функциональную зависимость между величинами x и y , которые получены в результате измерений.

При аналитическом исследовании взаимосвязи между двумя величинами x и y производят ряд наблюдений и в результате получается таблица значений:

x			¼		¼
y			¼		¼

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

Между величинами x и y существует функциональная зависимость, но ее аналитический вид обычно неизвестен, поэтому возникает практически важная задача - найти эмпирическую формулу

(1)

(где
- параметры), значения которой при
возможно мало отличались бы от опытных значений
.

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

Если в эмпирическую формулу (1) подставить исходные
, то получим теоретические значения
, где
.

Разности
называются отклонениями и представляют собой расстояния по вертикали от точек
до графика эмпирической функции.

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

(2)

будет минимальной.

Поясним геометрический смысл метода наименьших квадралтов.

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

Построение эмпирической формулы состоит из двух этапов: выяснение общего вида этой формулы и определение ее наилучших параметров.

Если неизвестен характер зависимости между данными величинами x и y, то вид эмпирической зависимости является произвольным. Предпочтение отдается простым формулам, обладающим хорошей точностью. Удачный выбор эмпирической формулы в значительной мере зависит от знаний исследователя в предметной области, используя которые он может указать класс функций из теоретических соображений. Большое значение имеет изображение полученных данных в декартовых или в специальных системах координат (полулогарифмической, логарифмической и т.д.). По положению точек можно примерно угадать общий вид зависимости путем установления сходства между построенным графиком и образцами известных кривых.

Определение наилучших коэффициентов , входящих в эмпирическую формулу производят хорошо известными аналитическими методами.

Для того, чтобы найти набор коэффициентов
, которые доставляют минимум функции S
, определяемой формулой (2), используем необходимое условие экстремума функции нескольких переменных - равенство нулю частных производных. В результате получим нормальную систему для определения коэффициентов
:

(3)

Таким образом, нахождение коэффициентов сводится к решению системы (3).

Эта система упрощается, если эмпирическая формула (1) линейна относительно параметров
, тогда система (3) - будет линейной.

Конкретный вид системы (3) зависит от того, из какого класса эмпирических формул мы ищем зависимость (1). В случае линейной зависимости
система (3) примет вид:

(4)

Эта линейная система может быть решена любым известным методом (методом Гаусса, простых итераций, формулами Крамера).

В случае квадратичной зависимости
система (.3) примет вид:

(5)

Метод Гаусса решения систем линейных алгебраических уравнений

Метод Гаусса состоит в последовательном исключении неизвестных до тех пор, пока не останется одно уравнение с одним неизвестным.

При этом матрица СЛАУ приводится треугольному виду, где ниже главной диагонали располагаются только нули.

Приведение матрицы к треугольному виду называется прямым ходом метода Гаусса. Обратный ход начинается с решения последнего уравнения и заканчивается определением первого неизвестного.

Имеем Ax=b (6), где A=[aij
] - матрица размерности n*n, det A0, b=(b1
, b2
, …, bn
)T
.

В предположении, что a11
0, первое уравнение системы (6):

делим на коэффициент a11
, в результате получаем уравнение

.

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

Первое неизвестное оказалось исключенным из всех уравнений, кроме первого. Далее предполагаем, что a1
22
0, делим второе уравнение на a1
22
и исключаем неизвестное x2
из всех уравнений, начиная со второго, и т.д. В результате последовательного исключения неизвестных система уравнений преобразуется в систему уравнений с треугольной матрицей:

(7)

Совокупность проведенных действий называется прямым ходом метода Гаусса.

Из n-го уравнения системы (6) определяем xn
, из (n-1)-го - xn-1
и т.д. до x1
. Совокупность таких действий называется обратным ходом метода Гаусса.

Реализация прямого хода требует арифметических операций, а обратного - арифметических операций.

Блок-схема алгоритма программы на языке Pascal

center;">Практическая часть

Текст программы на языке Pascal, перечень использованных в программе идентификаторов и полученные результаты

Листинг программы:

program
MultiChain;

const

xy: array
[1..40,1..2] of
real = ((0.15, 7.26), (-0.01, 7.95),

(0.46, 8.77), (1.71, 9.72), (3.94, 10.78), (7.40, 12.45), (12.41, 14.74), (19.47, 18.18), (29.21, 23.36), (42.49, 30.99), (60.43, 41.97), (84.42, 57.38), (116.18, 78.45), (157.80, 106.65), (211.78, 143.64), (281.04, 191.28), (369.00, 251.71), (479.59, 327.27), (617.30, 420.58), (787.21, 534.52), (995.04, 672.24), 1247.17, 837.2), (1550.71, 1033.14), (1913.5, 1264.13), (2344.18, 1534.54), (2852.21, 1849.11), (3447.91, 2212.91), (4142.52, 2631.36), (4948.19, 3110.26), (5878.09, 3665.8), (6946.36, 4274.55), (8768.24, 4973.5), (9560.03, 5760.63), (11139.19,6641.98), (12924.34,7627.6), (14935.29,8725.62), (17193.14,9945.22), (19720.24,11296.05), (22540.29,12788.24), (25678.32,14432.44));

var

matr: array
[1..20,1..20] of
real; { матрица решения }

sum: array
[1..21] of
real; { коэффициенты матрицы }

B: array
[1..20] of
real; { временный вектор }

A: array
[1..20] of
real;

n: integer;

m: integer;

procedure
FillMatr;

var

i,j,z: integer;

ss: real;

begin

for
i:=1 to
m+1 do
sum[i]:=0; { инициализируем переменные }

for
i:=1 to
20 do
B[i]:=0;

for
z:=1 to
40 do begin

ss:=1;

for
i:=1 to
m+1 do begin
{ подсчитываем коэффициенты - степени переменных }

sum[i] := sum[i] + ss;

ss := ss * xy[z,1];

end
;

ss:=xy[z,2];

for
i:=1 to
n+1 do begin

b[i] := b[i] + ss;

ss := ss * xy[z,1];

end
;

end
;

for
i:=1 to
n+1 do
matr[n+2,i] := b[i];

for
i:=1 to
n+1 do

for
j:=1 to
n+1 do

begin

matr[i,j] := sum[i + (j-1)];

end
;

end
;

procedure
TryFind;

var

i,j,z:integer;

Q: real;

begin

{1 step}

for
z:=1 to
n do begin

Q := matr[z,z];

for
j:=1 to
n+2 do

matr[j,z]:= matr[j,z] / Q;

for
i:=z+1 to
n+1 do begin

Q := (-1) / matr[z,i];

for
j:=1 to
n+2 do

matr[j,i]:= Q * matr[j,i] + matr[j,z];

end
;

end
;

Q := matr[n+1,n+1];

for
j:=n+1 to
n+2 do

matr[j,n+1]:= matr[j,n+1] / Q;

{2 step}

for
z:=n downto
1 do begin

Q := matr[n+2,z+1];

for
i:=z downto
1 do begin

matr[n+2,i]:=matr[n+2,i] - matr[z+1,i] * Q;

matr[z+1,i]:=0;

end
;

end
;

Writeln;

Write('F(x)=');

for
i:=1 to
n+1 do begin

A[i]:=matr[n+2,i];

if
i=1

then
Write('(',A[i]:0:4,')')

else
Write('(',A[i]:0:4,')*x^',i-1);

if
i<>n+1

then
Write('+')

else
Writeln

end
;

end
;

function
func(x:real):real;

{ функция для вычисления значения искомой функции по вычисленному полиному }

var

i: integer;

c,ff: real;

begin

c:=1;

ff:=0;

for
i:=1 to
n+1 do begin

ff := ff + A[i] * c;

c := c * x;

end
;

func := ff;

end
;

procedure
TryCheck;

var

i:Integer;

x,y,f,

delta:real;

begin

delta := 0;

for
i:=1 to
40 do begin

x:=xy[i,1]; y:=xy[i,2]; f:=func(x);

delta := delta + sqr(y - f);

end
;

delta := sqrt(1/40*delta);

Writeln('delta = ',delta:0:6);

end
;

begin

n := 2; { степень полинома }

writeln('Степень полинома n=',n);

m := 2*n; { количество коэффициентов }

FillMatr;

TryFind;

TryCheck;

writeln;

n := 4; { степень полинома }

writeln('Степень полинома n=',n);

m := 2*n; { количество коэффициентов }

FillMatr;

TryFind;

TryCheck;

writeln;

n := 5; { степень полинома }

writeln('Степень полинома n=',n);

m := 2*n; { количество коэффициентов }

FillMatr;

TryFind;

TryCheck;

end
.

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

xy: array [1..40,1..2] of real -
массив исходных значений;

matr: array [1..20,1..20] of real
- основная матрица для нахождения решения;

sum: array [1..21] of real -
массив коэффициентов для матрицы поиска;

B: array [1..20] of real -
массив вектора свободных членов;

A: array [1..20] of real -
массив для искомых коэффициентов полинома;

n: intege r -
переменная для размерности матрицы поиска;

m: integer -
переменная для размерность массива коэффициентов.

Переменные, встречающиеся в большинстве процедур:

i, j, z: integer
- переменные итераций для перебора элементов массива.

Переменные в процедурах:

Процедура FillMatr

ss: real
- переменная для подсчета коэффициентов матрицы;

Процедура TryFind

Q: real -
переменная для вычисления множителя;

Функция func

c: real -
переменная для вычисления степени переменной x(от 0-ой до n-ой);

ff: real
- переменная для вычисления значения функции-полинома;

Процедура TryCheck

x: real -
заданное значение независимой переменной;

y: real
- заданное значения неизвестной функции в точке x;

f: real -
переменная для значения вычисленной функции-полинома в точке x;

delta: real -
среднеквадратичное отклонение.

Программа, разработанная в среде PascalABC дает следующие результаты:

Описание решения задачи в среде MathCAD

Подберем полиномы второй, четвертой и пятой степени, в качестве аппроксимирующей функции. Для этих целей служат встроенные функции regress

и функция interp

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

Функция regress


является вспомогательной, она подготавливает дан­ные, необходимые для работы функции interp

. Вектор vs


содержит, в том числе, и коэффициенты полинома

Функция interp


возвращает значение полинома в точке z

. Определив новые функции f
2,
f
4
и f
5
, получим возможность находить значение поли­нома в любой заданной точке.

Коэффициенты задаем следующим образом:

Коэффициенты имеют вид:

Среднеквадратичные отклонения для каждого из трех случаев:

Результаты вычислений и их анализ

Представим результаты в виде сравнительных таблиц для применявшихся сред вычисления.

Для степени полинома n=2 получены следующие коэффициенты и среднеквадратичное отклонение (с точностью 0,001):

Среда вычисления	Коэффициенты			Среднеквадратичное отклонение
	А0	А1	А2
Pascal	32.266	0.618	0	61.597
MathCAD	32.264	0.618	0	61.597

Итак, коэффициенты полинома второй степени, полученные в двух средах вычислений фактически не отличаются (отличается лишь А0
на 0,002). Получен следующий полином второй степени, аппроксимирующий исходную табулируемую функцию:

.

Для степени полинома n=4 получены следующие коэффициенты и среднеквадратичное отклонение (с точностью 0,001):

Среда вычисления	Коэффициенты					Среднеквадратичное отклонение
	А0	А1	А2	А3	А4
Pascal	7,762	0,674	0	0	0	50,905
MathCAD	7,761	0,674	0	0	0	50,905

Итак, коэффициенты полинома четвертой степени, полученные в двух средах вычислений фактически не отличаются (отличается лишь А0
на 0,001). Получен следующий полином четвертой степени, аппроксимирующий исходную табулируемую функцию:

.

Для степени полинома n=5 получены следующие коэффициенты и среднеквадратичное отклонение (с точностью 0,001):

Среда вычисления	Коэффициенты						Среднеквадратичное отклонение
	А0	А1	А2	А3	А4	А5
Pascal	-0,213	0,71	0	0	0	0	47,9
MathCAD	-0,216	0,71	0	0	0	0	47,9

Итак, коэффициенты полинома четвертой степени, полученные в двух средах вычислений фактически не отличаются (отличается лишь А0
на 0,003). Получен следующий полином четвертой степени, аппроксимирующий исходную табулируемую функцию:

.

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

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

В среде MathCAD построены графики трех полученных приближающих функций в одной системе координат с исходными точками (х
i
,
yi
)
:

Рис.1. График 4-х функций в одной системе координат

Заключение

Аппроксимация функции методом наименьших квадратов является простой и легко реализуемой как в среде Pascal, так и в MathCAD. МНК «сглаживает» функцию, выбирая промежуточные значения, что является выгодным решением.

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

Результаты вычислений в среде Pascal и MathCAD отличаются не значительно. Однако выявлено, что при увеличении степени полинома, среднеквадратичное отклонение уменьшается, что означает, что построенный полином пятой степени ближе к точным исходным значениям табулированной функции.

Список литературы

1. А.А.Дадаян. Алгебра и геометрия./А.А.Дадаян, В.А.Дударенко. Минск: „Вышэйная школа”, 1989г

2. Вычислительная техника и программирование. Под ред. А.В. Петрова. М.: Высшая школа, 1991.

3. Информатика: Методические указания к курсовой работе. Санкт-Петербургский горный институт. Сост. Д.Е. Гусев, Г.Н. Журов. СПб, 1999

4. Индейкин В. В. Табличный редактор Microsoft Excel. Учебное пособие. – Казань, 1999. – 75с.

5. Кудрявцев Е. М. MathCAD 2000 Pro. – М.: ДМК Пресс, 2001. – 571с.

6. Крылов В.И., Бобков В.В., Монастырный П.И. Вычислительные методы:, т.2.-М.:Наука. Гл. ред. физ.-мат. лит., 1977.

7. Марьямов А. Г. "Применение модульного способа програмирования в среде Turbo Pascal 7.0 с целью решения полной задачи линейного программирования".

8. Моргун Александр Николаевич Программирование на языке Паскаль (Pascal). Основы обработки структур данных. — М.: «Диалектика», 2005. — С. 576. — ISBN 5-8459-0935-X.

9. Перминов Олег Николаевич Язык программирования Паскаль : Справочник. — М.: «Радио и связь», 1989. — С. 128. — ISBN 5-256-00311-9.