2. Аппроксимация в пакете MathCad

• Создание исходных массивов для работы

• Аппроксимация массива полиномами: 1, 5, 10, 15, 20, 22, 23 степени:

• Расчет среднеквадратичного отклонения для полиномов:

3. Реализация аппроксимации на яву:

• Принцип работы программы

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

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

Пусть на некотором отрезке в точках x₀, x₁, x₂… x_n нам известны значения y₀, y₁, y₂… y_n некоторой функции f(x).

Требуется определить параметры a_i многочлена вида F(x)=a₀+a₁x+a₂x₂+...+a_kx^k, где k<N

такого, что сумма квадратов отклонений значений y от значений функции f(y) в заданных точках x была минимальной, т.е.

S=Σ(y_i - f(x_i))² → min

Геометрический смысл заключается в том, что график найденного многочлена y = f(x) будет проходить как можно ближе к каждой из заданных точек.

Такая задача будет решена, если решить систему уравнений вида:

a₀n +a₁Σx_i+a₂Σx_i²+... +a_kΣx_i^k=Σy_i

a₀Σx_i +a₁Σx_i²+a₂Σx_i³+... +a_kΣx_i^k+1=Σx_iy_i

...

a₀Σx_i^k+a₁Σx_i^k+1+a₂Σx_i^k+2+...+a_kΣx_i^2k=Σx_i^ky_i

где везде под символом Σ подразумевается суммирование по i от 0 до n.

Для решения системы воспользуемся методом Гаусса (система уравнений приводится к треугольному виду, после чего последовательно вычисляются аn, … а0). Имея аналитический вид функции f, нетрудно найти её значение в любой точке и построить кривую.

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

uses crt,graph;

const razr=12;

ymin=50;

ymax=54;

type

mas=array[0..30] of real;

matr=array[0..30,0..30] of real;

var

gd,gm,code:integer;

i, j, k, n,i1,j1,di,dj: Integer;

Eur, Day: mas;

A,Ai:matr;

B,C:mas;

DM:matr;

input:text;

st:string;

y1,m:integer;

procedure ReadFile;

begin

assign

(input,'Vmlab.txt');

reset (input);

for i:=0 to 30 do readln (input,Day[i],Eur[i]);

close (input);

end;

procedure OutMas (M:mas;n:integer);

var i:integer;

begin

for i:=0 to n do writeln (M[i]:razr:3);

end;

procedure OutMatr (M:matr;n1,n2:integer);

var i,j:integer;

begin

for j:=0 to n2 do

begin

for i:=0 to n1 do write (M[i,j]:razr:3);

writeln;

end;

end;

function Step (x:real;y:integer):real;

var

st:real;

i:integer;

begin

st:=1;

for i:=1 to y do st:=st*x;

step:=st;

end;

function Det (M:matr;n:integer):real{extended};

var i,j:integer;

d:extended;

begin

if n=1 then d:=M[0,0];

if n=2 then d:=M[0,0]*M[1,1]-M[0,1]*M[1,0];

if n=3 then d:=(M[0,0]*M[1,1]*M[2,2])+(M[1,0]*M[2,1]*M[0,2])+(M[0,1]*M[1,2]*M[2,0])-(M[0,2]*M[1,1]*M[2,0])

-(M[0,1]*M[1,0]*M[2,2])-(M[0,0]*M[2,1]*M[1,2]);

Det:=D;

end;

function poly (x:real; m:byte):real;

var p:real;

i:integer;

begin

p:=0;

for i:=0 to m do p:=P+C[i]*Step(x,i);

poly:=P;

end;

begin

clrscr;

readfile;

{A}

writeln('Please,enter degree of polynom');

readln(m);

for i:=0 to m do

for j:=0 to m do

begin

A[i,j]:=0;

for k:=0 to 30 do A[i,j]:=A[i,j]+Step(Day[k],i+j);

end;

{}OutMatr(A,m,m);{}

{B}

for i:=0 to m do

begin

B[i]:=0;

for k:=0 to 30 do B[i]:=B[i]+Eur[k]*Step(Day[k],i);

end;

{}OutMas (B,m);{}

{A^-1}

if det(A,m+1)=0 then begin writeln ('Error: Det=0');

readkey;

halt;

end;

for i:=0 to m do

for j:=0 to m do

begin

{di:=0j:=0;}

for i1:=0 to m do

for j1:=0 to m do

begin

if i1>i then di:=1 else di:=0;

if

j1>j then dj:=1 else dj:=0;

if (i<>i1) and (j<>j1) then DM[i1-di,j1-dj]:=A[i1,j1];

end;

Ai[i,j]:=step(-1,i+j)*det(DM,m)/det(A,m+1);

end;

clrscr;

{}OutMatr(Ai,m,m);{}

{C}

for i:=0 to m do begin

C[i]:=0;

for j:=0 to m do C[i]:=C[i]+Ai[i,j]*B[j]

end;

{}OutMas (C,m);{}

clrscr;

gd:=detect;

initgraph(gd,gm,'');

setcolor(White);

setlinestyle(0,0,3);

line(40,420,620,420);

line(40,40,40,420);

setlinestyle(1,0,1);

if m=2 then

outtextxy(10,10,'degree=2')

else

outtextxy(10,10,'degree=1');

for i:=0 to 30 do

begin

line(40+19*round(Day[i]),420,40+19*round(Day[i]),60);

str(round(Day[i]),st);

outtextxy(37+19*round(Day[i]),430,st);

end;

y1:=30;

for i:=0 to 12 do

begin

y1:=y1+30;

line (40,y1,610,y1);

str((ymax-(ymax-ymin)*i/12):0:2,st);

outtextxy(1,y1-4,st);

end;

for i:=0 to 30 do

begin

pieslice(40+19*round(Day[i]),420-round((Eur[i]-Ymin)*360/(ymax-ymin)),0,360,3);

end;

setcolor (red);

setfillstyle(1,red);

for i:=40 to 610 do begin pieslice(i,420-round((poly((i-40)/19,m)-Ymin)*360/(ymax-ymin)),0,360,1);

end;

readln;

end.