Аппроксимация МНК (стр. 3 из 4)

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;

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;

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. Определив новые функции f2, f4 и f5, получим возможность находить значение полинома в любой заданной точке.

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

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

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

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

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

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

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

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

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

Среда вычисления	Коэффициенты					Среднеквадратичное отклонение
Среда вычисления	А₀	А₁	А₂	А₃	А₄	Среднеквадратичное отклонение
Pascal	7,762	0,674	0	0	0	50,905
MathCAD	7,761	0,674	0	0	0	50,905

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

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

Среда вычисления	Коэффициенты						Среднеквадратичное отклонение
Среда вычисления	А₀	А₁	А₂	А₃	А₄	А₅	Среднеквадратичное отклонение
Pascal	-0,213	0,71	0	0	0	0	47,9
MathCAD	-0,216	0,71	0	0	0	0	47,9

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

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

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