Аппроксимация функций. Выбор эмпирических формул

Шаг 3. В ячейку K2 вводим формулу =(A2-$F$49)^2.

Шаг 4. В ячейки K3:K26 эта формула копируется.

Шаг 5. В ячейку L2 вводим формулу =(B2-$F$50)^2.

Шаг 6. В ячейки L3:L26 эта формула копируется.

Шаг 7. В ячейку M2 вводим формулу =($D$37+$D$38*A2-B2)^2.

Шаг 8. В ячейки M3:M26 эта формула копируется.

Шаг 9. В ячейку N2 вводим формулу

=($I$38+$I$39*A2+$I$40*A2^2-B2)^2.

Шаг 10. В ячейки N3:N26 эта формула копируется.

Шаг 11. В ячейку O2 вводим формулу

=($G$47*EXP($G$46*A2)-B2)^2.

Шаг 12. В ячейки O3:O26 эта формула копируется.

 

Последующие шаги делаем с помощью автосуммирования .

 

Шаг 13. В ячейку J27 вводим формулу =СУММ(J2:J26).

Шаг 14. В ячейку K27 вводим формулу =СУММ(K2:K26).

Шаг 15. В ячейку L27 вводим формулу =СУММ(L2:L26).

Шаг 16. В ячейку M27 вводим формулу =СУММ(M2:M26).

Шаг 17. В ячейку N27 вводим формулу =СУММ(N2:N26).

Шаг 18. В ячейку O27 вводим формулу =СУММ(O2:O26).

Теперь проведем расчеты коэффициента корреляции по формуле

(только для линейной аппроксимации)

и коэффициента детерминированности по формуле . Результаты расчетов средствами Microsoft Excel представлены в таблице 8.

Таблица 8

Результаты расчета.

 

 

В таблице 8 в ячейке D53 записана формула =J27/(K27*L27)^(1/2).

В ячейке D54 записана формула =1- M27/L27.

В ячейке D55 записана формула =1- N27/L27.

В ячейке D56 записана формула =1- O27/L27.

Анализ результатов расчетов показывает, что квадратичная аппроксимация наилучшим образом описывает экспериментальные данные.

2.5. Построение графиков в Excel и использование функции ЛИНЕЙН.

Рассмотрим результаты эксперимента, приведенные в исследованном выше примере.

Исследуем характер зависимости в три  этапа:

• Построим график зависимости.
• Построим линию тренда ( ,  , ).
• Получим числовые характеристики коэффициентов этого уравнения.

 

 

 

 

 

 

 

 

Рис.4.1. График зависимости y от x 

Рис.4.2. График линейной аппроксимации

 

Рис.4.3. График квадратичной аппроксимации.

                   Рис.4.4. График экспоненциальной аппроксимации.

 

2.6.Программа на языке Pascal.

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

Текст программы:

program Kramer;

uses CRT;

const

n=25;

type

TArrayXY = array[1..2,1..n] of real;

TArray = array[1..n] of real;

var

SumX,SumY,SumX2,SumXY,SumX3,SumX4,SumX2Y,SumLnY,SumXLnY: real;

OPRlin,OPRkvadr,OPRa1,OPRa2,OPRa3:real;

a1lin,a2lin,a1kvadr,a2kvadr,a3kvadr,a1exp,a2exp,cexp:real;

Xsr,Ysr,S1,S2,S3,Slin,Skvadr,Sexp:real;

Kkor,KdetLin,KdetKvadr,KdetExp:real;

i:byte;

const

ArrayXY:TArrayXY=((12.85,12.32,11.43,10.59,10.21,9.65,9.63,9.22,8.44,8.07,7.74,7.32,7.08,6.87,5.23,5.02,4.65,4.53,3.24,2.55,1.86,1.76,1.11,0.99,0.72) , (154.77

145.59,108.37,100.76,98.32,81.43,80.97,79.04,61.76,60.54,55.86,47.63,48.03,36.85,25.65,24.98,22.87,20.32,9.06,6.23,3.91,3.22,1.22,1.10,0.53));

begin

ClrScr;

SumX:=0.0;

SumY:=0.0;

SumXY:=0.0;

SumX2:=0.0;

SumX3:=0.0;

SumX4:=0.0;

SumX2Y:=0.0;

SumLnY:=0.0;

SumXLnY:=0.0;

{ Вычисление сумм x, y, x*y, x^2, x^3, x^4, (x^2)*y, Ln(y), x*Ln(y) }

for i:=1 to n do

  begin

   SumX:=SumX+ArrayXY[1,i];

   SumY:=SumY+ArrayXY[2,i];

   SumXY:=SumXY+ArrayXY[1,i]*ArrayXY[2,i];

   SumX2:=SumX2+sqr(ArrayXY[1,i]);

   SumX3:=SumX3+ArrayXY[1,i]*ArrayXY[1,i]*ArrayXY[1,i];

   SumX4:=SumX4+sqr(ArrayXY[1,i])*sqr(ArrayXY[1,i]);

   SumX2Y:=SumX2Y+sqr(ArrayXY[1,i])*ArrayXY[2,i];

   SumLnY:=SumLnY+ln(ArrayXY[2,i]);

   SumXLnY:=SumXLnY+ArrayXY[1,i]*ln(ArrayXY[2,i])

  end;

{ Вычисление коэффициентов }

 OPRlin:=0.0;

a1lin:=0.0;

a2lin:=0.0;

a1kvadr:=0.0;

OPRkvadr:=0.0;

a2kvadr:=0.0;

a2kvadr:=0.0;

a1exp:=0.0;

a2exp:=0.0;

OPRlin:=n*SumX2-SumX*SumX;

a1lin:=(SumX2*SumY-SumX*SumXY)/OPRlin;

a2lin:=(n*SumXY-SumX*SumY)/OPRlin;

OPRkvadr:=n*SumX2*SumX4+SumX*SumX3*SumX2+SumX2*SumX*SumX3-     SumX2*SumX2*SumX2-n*SumX3*SumX3-SumX*SumX*SumX4;

a1kvadr:=(SumY*SumX2*SumX4+SumX*SumX2Y*SumX3+SumX2*SumXY*SumX3- SumX2*SumX2*SumX2Y-SumY*SumX3*SumX3-SumX*SumXY*SumX4)/OPRkvadr;

a2kvadr:=(n*SumXY*SumX4+SumY*SumX3*SumX2+SumX2*SumX*SumX2Y-SumX2*SumX2*SumXY-n*SumX3*SumX2Y-SumY*SumX*SumX4)/OPRkvadr;

a3kvadr:=(n*SumX2*SumX2Y+SumX*SumXY*SumX2+SumY*SumX*SumX3-SumY*SumX2*SumX2-n*SumXY*SumX3-SumX*SumX*SumX2Y)/OPrkvadr;

a2exp:=(n*SumXLnY-SumX*SumLnY)/OPRlin;

cexp:=(SumX2*SumLnY-SumX*SumXLnY)/OPRlin;

a1exp:=exp(cexp);

 { Вычисление средних арифметических x и y }

 Xsr:=SumX/n;

Ysr:=SumY/n;

S1:=0.0;

S2:=0.0;

S3:=0.0;

Slin:=0.0;

Skvadr:=0.0;

Sexp:=0.0;

Kkor:=0.0;

KdetLin:=0.0;

KdetKvadr:=0.0;

KdetExp:=0.0;

 

for i:=1 to n do

  begin

   S1:=S1+(ArrayXY[1,i]-Xsr)*(ArrayXY[2,i]-Ysr);

   S2:=S2+sqr(ArrayXY[1,i]-Xsr);

   S3:=S3+sqr(ArrayXY[2,i]-Ysr);

   Slin:=Slin+sqr(a1lin+a2lin*ArrayXY[1,i]-ArrayXY[2,i]);

Skvadr:=Skvadr+sqr(a1kvadr+a2kvadr*ArrayXY[1,i]+a3kvadr*ArrayXY[1,i]*ArrayXY[1,i]-ArrayXY[2,i]);

   Sexp:=Sexp+sqr(a1exp*exp(a2exp*ArrayXY[1,i])-ArrayXY[2,i]);

  end;

{ Вычисление коэффициентов корреляции  и детерминированности }

 Kkor:=S1/sqrt(S2*S3);

KdetLin:=1-Slin/S3;

KdetKvadr:=1-Skvadr/S3;

KdetExp:=1-Sexp/S3;

 { Вывод результатов }

 WriteLn('Линейная функция');

 WriteLn('a1=',a1lin:8:5);

WriteLn('a2=',a2lin:8:5);

 WriteLn('Квадратичная функция');

 WriteLn('a1=',a1kvadr:8:5);

WriteLn('a2=',a2kvadr:8:5);

WriteLn('a3=',a3kvadr:8:5);

 WriteLn('Экспоненциальная функция');

 WriteLn('a1=',a1exp:8:5);

WriteLn('a2=',a2exp:8:5);

WriteLn('c=',cexp:8:5);

WriteLn('Xcp=',Xsr:8:5);

WriteLn('Ycp=',Ysr:8:5);

WriteLn('Коэффициент корреляции ',Kkor:8:5);

 WriteLn('Коэффициент детерминированности (линейная аппроксимация) ',KdetLin:2:5);

WriteLn('Коэффициент детерминированности (квадратическая аппроксимация) ',KdetKvadr:2:5);

WriteLn('Коэффициент детерминированности (экспоненциальная аппроксимация) ',KdetExp:2:5);

end.

Результаты расчета Pascal.

Коэффициенты линейной функции

a1=-24.73516

a2=11.63471

Коэффициенты квадратичной функции

a1= 1.59678

a2=-0.62145

a3= 0.95543

Коэффициенты экспоненциальной функции

a1= 1.65885

a2= 0.40987

c= 0.50613

Xcp= 6.52320

Ycp=51.16040

Коэффициент корреляции  0.96196

Коэффициент детерминированности (линейная аппроксимация) 0.92537

Коэффициент детерминированности (квадратическая аппроксимация) 0.99409

Коэффициент детерминированности (экспоненциальная аппроксимация) 0.02691

 

 

 

 

 

Схема алгоритма.

Рис.5.1. Блок-схема

 

 

 Результаты выполнения программы:

 

Заключение.

Сделаем заключение по результатам полученных данных:

1. Анализ  результатов расчетов показывает, что квадратичная аппроксимация наилучшим образом описывает экспериментальные данные т.к. согласно таблице 8 коэффициент корреляции - 0,9620; Коэффициенты детерминированности линейной аппроксимации - 0,9253; квадратической аппроксимации – 0,994; экспоненциальной аппроксимация – 0,0269.

2. Сравнивая  результаты, полученные при помощи  функции ЛИНЕЙН видим что они  полностью совпадают с вычислениями, проведенными выше. Это указывает  на то, что вычисления верны.

3. Полученное  при построении линии тренда  значение коэффициента детерминированности для экспоненциальной зависимости не совпадает с истинным значением поскольку при вычислении коэффициента детерминированности используются не истинные значения y, а преобразованные значения ln(y) с дальнейшей линеаризацией.

4. Результаты  полученные с помощью программы на языке PASCAL и Delphi полностью совпадают со значениями приведенными выше. Это говорит о верности вычислений.

 

 

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

1. Ахметов К.С. Windows 95 для всех. - М.:ТОО "КомпьютерПресс", 1995.
2. Вычислительная техника и программирование. Под ред. А.В. Петрова. М.: Высшая школа, 1991.             
3. Гончаров A., Excel 97 в примерах. — СПб: Питер, 1997.
4. Левин А., Самоучитель работы на компьютере. - М.: Международное агентство А.Д.Т., 1996.
5. Информатика: Методические указания к курсовой работе. Санкт-Петербургский горный институт. Сост. Д.Е. Гусев, Г.Н. Журов. СПб, 1999

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

ПРИЛОЖЕНИЕ

unit Unit1;

 

interface

 

uses

  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,

  Dialogs, Math, StdCtrls, Grids, Menus, ComCtrls, Buttons, TeEngine,

  Series, ExtCtrls, TeeProcs, Chart, Unit2;

 

type

  TForm1 = class(TForm)

    Label1: TLabel;

    Label2: TLabel;

    Label3: TLabel;

    Label4: TLabel;

    Label5: TLabel;

    Label6: TLabel;

    Label7: TLabel;

    Label8: TLabel;

    Label9: TLabel;

    Label10: TLabel;

    Label11: TLabel;

    Label12: TLabel;

    Label13: TLabel;

    Edit1: TEdit;

    Edit2: TEdit;

    Edit3: TEdit;

    Edit4: TEdit;

    Edit5: TEdit;

    Edit6: TEdit;

    Edit7: TEdit;

    Edit8: TEdit;

    Edit9: TEdit;

    Edit10: TEdit;

    Edit11: TEdit;

    Edit12: TEdit;

    Edit13: TEdit;

    StatusBar1: TStatusBar;

    MainMenu1: TMainMenu;

    N1: TMenuItem;

    StringGrid1: TStringGrid;

    Label14: TLabel;

    Label15: TLabel;

    Label16: TLabel;

    Label17: TLabel;

    Label18: TLabel;

    Label19: TLabel;

    Label20: TLabel;

    Label21: TLabel;

    Label22: TLabel;

    Edit14: TEdit;

    Edit15: TEdit;

    Edit16: TEdit;

    Edit17: TEdit;

    Edit18: TEdit;

    Edit19: TEdit;

    Edit20: TEdit;

    Edit21: TEdit;

    Edit22: TEdit;

    BitBtn1: TBitBtn;

    Label23: TLabel;

    Label24: TLabel;

    Edit23: TEdit;

    Edit24: TEdit;

    Chart1: TChart;

    Series1: TLineSeries;

    Series2: TLineSeries;

    N2: TMenuItem;

    N3: TMenuItem;

    N4: TMenuItem;

    N5: TMenuItem;

    N6: TMenuItem;

    GroupBox1: TGroupBox;

    Label25: TLabel;

    Label26: TLabel;

    procedure FormCreate(Sender: TObject);

    procedure BitBtn1Click(Sender: TObject);

    procedure N4Click(Sender: TObject);

    procedure N6Click(Sender: TObject);

  private

    { Private declarations }

  public

    { Public declarations }

  end;

  mas=array[1..9]of real;

 

var

  Form1: TForm1;

  n:byte;

  x,y:mas;

 

implementation

 

//--Раздел условия--------------------------------------------------------------

procedure sets(var x,y:mas);

begin

//Функция, заданная таблично

x[1]:=1;   y[1]:=1.24;

x[2]:=1.5; y[2]:=1.74;

x[3]:=2;   y[3]:=1.61;

x[4]:=2.5; y[4]:=2.16;

x[5]:=3;   y[5]:=3.06;

x[6]:=3.5; y[6]:=2.88;

x[7]:=4;   y[7]:=4.53;

x[8]:=4.5; y[8]:=5.4;

x[9]:=5;   y[9]:=7.07;

end;

//------------------------------------------------------------------------------

 

procedure approximation(x,y:mas; n:byte; var xa,xh,xg,ya,yh,yg,zxa,zxh,zxg:real; var delta:mas; var deltak,a,b,kvadra:real; var zz:string; var zzz:mas; var aa0,bb0,cc0,dd10,dd20:real);

var p:real;

    k,i:byte;

    aaa,bbb,aa,bb,cc,dd1,dd2:mas;

 

function z(k:byte; a,b,x:real):real;

begin

case k of

1: Result:=a*x+b;

2: Result:=b*Power(x,a);

3: Result:=b*exp(a*x);

4: Result:=a*ln(x)+b;

5: Result:=a/x+b;

6: Result:=1/(a*x+b);

7: Result:=x/(b*x+a);

8: Result:=b*exp(a/x);

9: Result:=1/(a*ln(x)+b);

end;

end;

function fi1(k:byte; x:real):real;

begin

case k of

1: Result:=x;

2: Result:=ln(x);

3: Result:=x;

4: Result:=ln(x);

5: Result:=1/x;

6: Result:=x;

7: Result:=1/x;

8: Result:=1/x;

9: Result:=ln(x);

end;

end;

function zx(k:byte; xa,xh,xg:real):real;

begin

case k of

1: Result:=xa;

2: Result:=xg;

3: Result:=xa;

4: Result:=xg;

5: Result:=xh;

6: Result:=xa;

7: Result:=xh;

8: Result:=xh;

9: Result:=xg;

end;

end;

function dd(k,i:byte; y:mas):real;

begin

case k of

1: Result:=y[i];

2: Result:=ln(y[i]);

3: Result:=ln(y[i]);

4: Result:=y[i];

5: Result:=y[i];

6: Result:=1/y[i];

7: Result:=1/y[i];

8: Result:=ln(y[i]);

9: Result:=1/y[i];

end;

end;

function zzx(k:byte; a,b:real):string;

begin

case k of

1: Result:=FloatToStrF(a,fffixed,6,2)+'*x+'+FloatToStrF(b,fffixed,6,2);

2: Result:=FloatToStrF(b,fffixed,6,2)+'*Power(x,'+FloatToStrF(a,fffixed,6,2)+')';

3: Result:=FloatToStrF(b,fffixed,6,2)+'*exp('+FloatToStrF(a,fffixed,6,2)+'*x)';

4: Result:=FloatToStrF(a,fffixed,6,2)+'*ln(x)+'+FloatToStrF(a,fffixed,6,2);

5: Result:=FloatToStrF(a,fffixed,6,2)+'/x+'+FloatToStrF(b,fffixed,6,2);

6: Result:='1/('+FloatToStrF(a,fffixed,6,2)+'*x+'+FloatToStrF(b,fffixed,6,2)+')';

7: Result:='x/('+FloatToStrF(b,fffixed,6,2)+'*x+'+FloatToStrF(a,fffixed,6,2)+')';

8: Result:=FloatToStrF(b,fffixed,6,2)+'*exp('+FloatToStrF(a,fffixed,6,2)+'/x)';

9: Result:='1/('+FloatToStrF(a,fffixed,6,2)+'*ln(x)+'+FloatToStrF(a,fffixed,6,2)+')';

end;

end;

 

begin

xa:=(x[1]+x[n])/2;

xh:=sqrt(x[1]*x[n]);

xg:=(2*x[1]*x[n])/(x[1]+x[n]);

ya:=(y[1]+y[n])/2;

yh:=sqrt(y[1]*y[n]);

yg:=(2*y[1]*y[n])/(y[1]+y[n]);

//Находим все delta 1..9

for k:=1 to 9 do

  begin

  aa[k]:=0;

  for i:=1 to n do

    aa[k]:=aa[k]+(fi1(k,x[i])*fi1(k,x[i]));

  bb[k]:=0;

  for i:=1 to n do

    bb[k]:=bb[k]+fi1(k,x[i]);

  cc[k]:=n;

  dd1[k]:=0;

  for i:=1 to n do

    dd1[k]:=dd1[k]+(fi1(k,x[i])*dd(k,i,y));

  dd2[k]:=0;

  for i:=1 to n do

    dd2[k]:=dd2[k]+dd(k,i,y);

  bbb[k]:=((aa[k]*dd2[k])-(bb[k]*dd1[k]))/((aa[k]*cc[k])-(bb[k]*bb[k]));

  aaa[k]:=(dd1[k]-bb[k]*bbb[k])/aa[k];

  delta[k]:=Abs(z(k,aaa[k],bbb[k],zx(k,xa,xh,xg)));

  end;

//Находим deltak

deltak:=delta[1];

k:=1;

for i:=2 to 9 do

  if delta[i]<deltak then

    begin

    deltak:=delta[i];

    k:=i;

    end;

a:=aaa[k];

if k in [2,3,8] then

  b:=exp(bbb[k])

  else

  b:=bbb[k];

aa0:=aa[k];

bb0:=bb[k];

cc0:=cc[k];

dd10:=dd1[k];

dd20:=dd2[k];

zz:=zzx(k,a,b);

for i:=1 to n do

  zzz[i]:=z(k,a,b,x[i]);

p:=0;

for i:=1 to n do

  p:=p+((zzz[i]-y[i])*(zzz[i]-y[i]));

kvadra:=sqrt(p);

zxa:=z(k,a,b,xa);

zxh:=z(k,a,b,xh);

zxg:=z(k,a,b,xg);

end;

 

{$R *.dfm}

 

procedure TForm1.FormCreate(Sender: TObject);

var i:byte;

begin

n:=StrToInt(Edit1.Text);

sets(x,y);

StringGrid1.RowCount:=n+1;

StringGrid1.Cells[1,0]:='x';

StringGrid1.Cells[2,0]:='y';

StringGrid1.Cells[3,0]:='z';

for i:=1 to n do

  begin

  StringGrid1.Cells[0,i]:=IntToStr(i);

  StringGrid1.Cells[1,i]:=FloatToStrF(x[i],fffixed,6,4);

  StringGrid1.Cells[2,i]:=FloatToStrF(y[i],fffixed,6,4);

  StringGrid1.Cells[3,i]:='';

  end;

Series1.Clear;

Series2.Clear;

for i:=1 to n do

  Series1.AddXY(x[i],y[i]);

end;

 

procedure TForm1.BitBtn1Click(Sender: TObject);

var xa,xh,xg,ya,yh,yg,zxa,zxh,zxg,kvadra,deltak,a,b,aa,bb,cc,dd1,dd2:real;

    z:string;

    delta,zz:mas;

    i:byte;

begin

approximation(x,y,n,xa,xh,xg,ya,yh,yg,zxa,zxh,zxg,delta,deltak,a,b,kvadra,z,zz,aa,bb,cc,dd1,dd2);

for i:=1 to n do

  StringGrid1.Cells[3,i]:=FloatToStrF(zz[i],fffixed,6,4);

Edit2.Text:=FloatToStrF(xa,fffixed,6,2);

Edit3.Text:=FloatToStrF(xh,fffixed,6,2);

Edit4.Text:=FloatToStrF(xg,fffixed,6,2);

Edit5.Text:=FloatToStrF(ya,fffixed,6,2);

Edit6.Text:=FloatToStrF(yh,fffixed,6,2);

Edit7.Text:=FloatToStrF(yg,fffixed,6,2);

Edit8.Text:=FloatToStrF(zxa,fffixed,6,2);