Algoritmy numerickej aproximácie (riešené príklady): Rozdiel medzi revíziami

Z Kiwiki
Skočit na navigaci Skočit na vyhledávání
Riadok 336: Riadok 336:
 
===Riešenie v jazyku C===
 
===Riešenie v jazyku C===
 
Kriteriom na vyhodnotenie najlepšej aproximovanej krivky bude hľadanie minimálych hodnôt a vzdialeností uzlových (vstupných) bodov od aproximujúcej krivky.
 
Kriteriom na vyhodnotenie najlepšej aproximovanej krivky bude hľadanie minimálych hodnôt a vzdialeností uzlových (vstupných) bodov od aproximujúcej krivky.
 +
 +
Definovanie štruktúr potrebných pre spracovanie súradníc:
 +
<source lang="c" line>
 +
struct Bod{
 +
  double x,y;
 +
};
 +
 +
enum TYP_KRIVKY{
 +
  LIN, LOG, EXP, MOC
 +
  };
 +
 +
struct Krivka{
 +
  TYP_KRIVKY typ;
 +
  double a,b;
 +
};
 +
</source>
 +
 +
Táto štruktúra pre zjednodušenie poslúži na uloženie všetkých typov aproximácie do jednej premennej typu Aprox:
 +
<source lang="c" line>
 +
struct Aprox{
 +
    Krivka Lin, Log, Exp, Moc;
 +
};
 +
</source>
 +
 +
Funkcia aproximacie vypočíta všetky 4 typy aproximácií. Návratová hodnota funkcie je typu Aprox:
 +
<source lang="c" line>
 +
Aprox aproximacie(Bod *body, int n, Aprox B)  //vypocita vsetky
 +
                                                //aproximacie
 +
{
 +
  double sx=0,sxx=0,sy=0,sxy=0,slxly;
 +
  double slx=0,slxx=0,sly=0,slxy=0,sxly=0;
 +
for(int i=0;i<n;i++)
 +
{    sx+=body[i].x;
 +
slx+=log(body[i].x);
 +
sy+=body[i].y;
 +
sly+=log(body[i].y);
 +
sxy+=body[i].x * body[i].y;
 +
slxy+=log(body[i].x) * body[i].y;
 +
sxly+=body[i].x * log(body[i].y);
 +
slxly+=log(body[i].x) *log(body[i].y);
 +
        sxx+=body[i].x * body[i].x;
 +
        slxx+=log(body[i].x) * log(body[i].x);
 +
    }
 +
 +
B.Lin.a=(n*sxy-sy*sx)/(n*sxx-sx*sx);
 +
B.Lin.b=(sy-B.Lin.a*sx)/n;
 +
 +
B.Log.a=(n*slxy-sy*slx)/(n*slxx-slx*slx);
 +
B.Log.b=(sy-B.Log.a*slx)/n;
 +
 +
B.Exp.a=(n*sxly-sly*sx)/(n*sxx-sx*sx);
 +
B.Exp.b=(sly-B.Exp.a*sx)/n;
 +
B.Exp.b=exp(B.Exp.b);
 +
 +
B.Moc.a=(n*slxly-sly*slx)/(n*slxx-slx*slx);
 +
B.Moc.b=(sly-B.Moc.a*slx)/n;
 +
    B.Moc.b=exp(B.Moc.b);
 +
 +
    return B;
 +
}
 +
</source>
 +
 +
Na vypísanie najvhodnejšieho typu rovnice je použitá funkcia vypisRovnicuKrivky:
 +
<source lang="c" line>
 +
void vypisRovnicuKrivky(Krivka k)    //vypise rovnicu urciteho typu
 +
{
 +
  cout<<"y = ";
 +
  switch(k.typ)
 +
  {
 +
  case LIN:
 +
  cout<<k.a<<" x";
 +
  if(k.b>0) cout<<"+";
 +
  cout<<k.b<<"  -> LINEARNA"<<endl;
 +
  break;
 +
  case LOG:
 +
  cout<<k.a<<" ln x";
 +
  if(k.b>0) cout<<"+";
 +
  cout<<k.b<<"  -> LOGARITMICKA"<<endl;
 +
  break;
 +
  case EXP:
 +
  cout<<k.b<<" e^(";
 +
  cout<<k.a<<" x)"<<"  -> EXPONENCIALNA"<<endl;
 +
  break;
 +
  case MOC:
 +
  cout<<k.b<<" x^(";
 +
  cout<<k.a<<" )"<<"  -> MOCNINOVA"<<endl;
 +
  break;
 +
}
 +
}
 +
</source>
 +
 +
Pre zistenie optimálnej aproximovanej krivky pre dané vstupné body je použitý jednoduchý algoritmus daný vo funkcii zistiOptimum.
 +
*Postup riešenia algoritmu spočíva v zistení vzdialenosti y-onovej súradnice vstupného bodu od aproximovanej y-onovej súradnice, ktorú musíme najskôr vypočítať, a porovnaním ako sú tieto body vzdialené od jednotlivých bodov kriviky pre každý typ, zistíme optimum. 
 +
*Dosadením x-ových súradníc vstupných bodov do príslušnej rovnice krivky sa vypočíta aproximovaná y-onová súradnica.  Absolútnu hodnotu rozdielu vzdialenosti medzi y-onovými súradnicami(vypočítanými a vstupnými) sčíta a uloží túto sumu do premennej sumLin (pre každý ďalší jeden typ aproximácie do - sumLog, sumExp, sumMoc ). Toto sa opakuje pre každý jeden vstupný bod a pre každú rovnicu krivky.
 +
*Na záver sa sumy všetkých typov aproximácií porovnajú. Najmenšia suma predstavuje najmenšie odchýlky pre danú krivku v porovnaní so vstupnými bodmi. Funkcia vráti Optimálnu krivku  aproximácie s najmenšou sumou:
 +
<source lang="c" line>
 +
Krivka zistiOptimum(Bod *body, int n, Aprox B)
 +
{
 +
double sumLin=0, sumLog=0, sumExp=0, sumMoc=0, y1=0,y2=0,y3=0,y4=0, Naj=0;
 +
Krivka Optim;
 +
for(int i=0; i < n; i++)
 +
{
 +
y1=((B.Lin.a*body[i].x)+B.Lin.b)-body[i].y; 
 +
if(y1<0) y1=-y1;
 +
  sumLin+=y1;
 +
 +
y2=((B.Log.a*(log(body[i].x)))+B.Log.b)-body[i].y; 
 +
                if(y2<0) y2=-y2;                         
 +
  sumLog+=y2;
 +
 +
y3=(B.Exp.b*(exp(B.Exp.a*body[i].x)));   
 +
y3=log(y3);
 +
if(y3<0) y3=-y3;
 +
sumExp+=y3;
 +
 +
y4=(B.Moc.b*(pow(body[i].x, B.Moc.a)))-body[i].y;
 +
if(y4<0) y4=-y4;                         
 +
  sumMoc+=y4;
 +
 +
}
 +
Naj=min(sumLin,sumLog);
 +
Naj=min(Naj,sumExp);
 +
Naj=min(Naj,sumMoc);
 +
 +
if(Naj==sumLin) {B.Lin.typ=LIN; return B.Lin;}
 +
if(Naj==sumLog) {B.Log.typ=LOG; return B.Log;}
 +
if(Naj==sumExp) {B.Exp.typ=EXP; return B.Exp;}
 +
 +
B.Moc.typ=MOC;
 +
return B.Moc;
 +
}
 +
</source>
 +
 +
===Modifikácia zdrojového kódu===
 +
Pre názornosť nebudeme v programe dáta načítavať z klávesnice, ale vstupné údaje budú zadané priamo v zdrojovom kóde
 +
<source lang="c" line>
 +
#include <iostream.h>
 +
#include <conio.h>
 +
#include <math.h>
 +
#include <cmath>
 +
 +
struct Bod{
 +
  double x,y;
 +
};
 +
enum TYP_KRIVKY{
 +
  LIN, LOG, EXP, MOC
 +
  };
 +
struct Krivka{
 +
  TYP_KRIVKY typ;
 +
  double a,b;
 +
};
 +
struct Aprox{
 +
Krivka Lin, Log, Exp, Moc;
 +
};
 +
 +
//Aprox aproximacie(Bod *body, int n, Aprox B)
 +
//void vypisRovnicuKrivky(Krivka k)
 +
//Krivka zistiOptimum(Bod *body, int n, Aprox B)
 +
 +
int main(int argc, char* argv[])
 +
{
 +
const int n=8;
 +
Bod body[n]={
 +
  {1.50,0.58},
 +
  {2.00,8.61},
 +
  {3.00,60.91},
 +
  {4.00,201.40},
 +
  {5.00,476.59},
 +
  {6.00,892.57},
 +
  {7.00,1304.37},
 +
  {8.00,2715.04}
 +
};
 +
Aprox B;
 +
cout<<"Najvhodnejsia: ";
 +
vypisRovnicuKrivky(zistiOptimum(body, n, B));
 +
 +
    getch();
 +
    return 0;
 +
}
 +
</source>

Verzia zo dňa a času 11:18, 27. máj 2010

Imbox draft.png
Toto je projekt, na ktorom sa ešte stále pracuje!!

Aj keď sú v tomto dokumente použiteľné informácie, ešte nie je dokončený. Svoje návrhy môžete vyjadriť v diskusii o tejto stránke.

Algoritmy a programovanie - zbierka úloh


Štruktúry

Rekurzia

Dynamická alokácia pamäti

Vyhľadávanie

Triedenie

Lineárny zoznam

Binárny strom

Numerické algoritmy

>Algoritmy numerickej interpolácie
>Algoritmy numerickej aproximácie
>Numerické integrovanie
>Numerické derivovanie

Zadanie 1

Riešte problém aproximácie dát. K dispozícii máme n bodov v rovine (ich súradnice x a y). Úlohou bude vypočítať rovnicu aproximujúcej krivky metódou najmenších štvorcov pre:

  1. Lineárnu aproximáciu v tvare [math]y=ax+b[/math]
  2. Logaritmickú aproximáciu v tvare [math]y=a \ln{x} + b[/math]
  3. Exponenciálnu aproximáciu v tvare [math]y=b e^{ax}[/math]
  4. Mocninovú aproximáciu v tvare [math]y=bx^a[/math]

Porovnajte výsledok funkcie s výsledkom v tabuľkovom procesore, kde do (bodového) grafu zo vstupných bodov pridáte spojnicu trendu typu lineárna, resp. logaritmická, mocninová exponenciálna.

Vstupné údaje

Ako vzorku vstupných bodov budeme uvažovať body z intervalu <1,8>. Pre generovanie vstupných hodnôt použijem funkciu [math]f(x)=x^4-exp(x)[/math]

Vzorový vstup

[math]x_i[/math] [math]y_i[/math]
1.50 0.58
2.00 8.61
3.00 60.91
4.00 201.40
5.00 476.59
6.00 892.57
7.00 1304.37
8.00 1115.04

Vzorový výstup

  1. y = 213.634 x - 467.195
  2. y = 784.982 ln x - 572.83
  3. y = 0.9628 e^(1.054x)
  4. y = 0.2609 x^(4.443)

Analýza matematických vzťahov

Poznámka: v tomto texte nie sú uvádzané celé postupy výpočtu, nakoľko v tomto predmete nie sú primárnym cieľom matematické dôkazy. Študent si postup môže ľahko overiť sám.

Lineárna aproximácia: y=ax+b

Pre lineárnu aproximáciu platia vzťahy:

[math]\begin{align} & b\sum\limits_{i=1}^{n}{{{x}_{i}}}+a\sum\limits_{i=1}^{n}{{{x}^{2}}_{i}}=\sum\limits_{i=1}^{n}{{{x}_{i}}{{y}_{i}}} \\ & nb+a\sum\limits_{i=1}^{n}{{{x}_{i}}}=\sum\limits_{i=1}^{n}{{{y}_{i}}} \\ \end{align}[/math]

Čo môžeme pomocou maticového zápisu zapísať nasledovne:

[math]\left[ \begin{matrix} \sum\limits_{i=1}^{n}{{{x}^{2}}_{i}} & \sum\limits_{i=1}^{n}{{{x}_{i}}} \\ \sum\limits_{i=1}^{n}{{{x}_{i}}} & n \\ \end{matrix} \right]\cdot \left[ \begin{matrix} a \\ b \\ \end{matrix} \right]=\left[ \begin{matrix} \sum\limits_{i=1}^{n}{{{x}_{i}}{{y}_{i}}} \\ \sum\limits_{i=1}^{n}{{{y}_{i}}} \\ \end{matrix} \right][/math]

Riešenie tejto sústavy rovníc je:

[math]\begin{align} & a=\frac{n\sum\limits_{i=1}^{n}{{{x}_{i}}{{y}_{i}}}-\sum\limits_{i=1}^{n}{{{x}_{i}}}\sum\limits_{i=1}^{n}{{{y}_{i}}}}{n\sum\limits_{i=1}^{n}{x_{i}^{2}}-{{\left( \sum\limits_{i=1}^{n}{{{x}_{i}}} \right)}^{2}}} \\ & b=\frac{\sum\limits_{i=1}^{n}{{{y}_{i}}}-a\sum\limits_{i=1}^{n}{{{x}_{i}}}}{n} \\ \end{align}[/math]

Logaritmická aproximácia: [math]y=a*\ln x + b[/math]

Hľadáme najlepšiu aproximáciu pre logaritmickú rovnicu: [math]y=a*\ln x + b[/math]. Postup výpočtu je rovnaký ako v prvom prípade. Jediná zmena v riešení je tá, že os x je v logaritmickej mierke. Dostávame teda sústavu 2 rovníc:

[math]\left[ \begin{matrix} \sum\limits_{i=1}^{n}{{{\left( \ln {{x}_{i}} \right)}^{2}}} & \sum\limits_{i=1}^{n}{\left( \ln {{x}_{i}} \right)} \\ \sum\limits_{i=1}^{n}{\left( \ln {{x}_{i}} \right)} & n \\ \end{matrix} \right]\cdot \left[ \begin{matrix} a \\ b \\ \end{matrix} \right]=\left[ \begin{matrix} \sum\limits_{i=1}^{n}{\left( \ln {{x}_{i}} \right){{y}_{i}}} \\ \sum\limits_{i=1}^{n}{{{y}_{i}}} \\ \end{matrix} \right][/math]

Riešenie tejto sústavy rovníc je:

[math]\begin{align} & a=\frac{n\sum\limits_{i=1}^{n}{\left( \ln {{x}_{i}} \right){{y}_{i}}}-\sum\limits_{i=1}^{n}{\ln {{x}_{i}}}\sum\limits_{i=1}^{n}{{{y}_{i}}}}{n\sum\limits_{i=1}^{n}{{{\left( \ln {{x}_{i}} \right)}^{2}}}-{{\left( \sum\limits_{i=1}^{n}{\ln {{x}_{i}}} \right)}^{2}}} \\ & b=\frac{\sum\limits_{i=1}^{n}{{{y}_{i}}}-a\sum\limits_{i=1}^{n}{\ln {{x}_{i}}}}{n} \\ \end{align}[/math]

Exponenciálna aproximácia: [math]y=b{{e}^{ax}}[/math]

Rovnicu [math]y=b{{e}^{ax}}[/math] si upravíme tak, že ju zapíšeme v logaritmickom tvare:

[math]\begin{align} & \ln y=\ln \left( b{{e}^{ax}} \right) \\ & \ln y=\ln b+\ln {{e}^{ax}} \\ & \ln y=ax+\ln b\,\,\left[ B=\ln b \right] \\ & \ln y=ax+B \\ \end{align}[/math]

Opäť dostávame sústavu 2 rovníc o 2 neznámych. V rovniciach je použitá substitúcia [math]B=\ln b[/math].

[math]\left[ \begin{matrix} \sum\limits_{i=1}^{n}{{{x}^{2}}_{i}} & \sum\limits_{i=1}^{n}{{{x}_{i}}} \\ \sum\limits_{i=1}^{n}{{{x}_{i}}} & n \\ \end{matrix} \right]\cdot \left[ \begin{matrix} a \\ B \\ \end{matrix} \right]=\left[ \begin{matrix} \sum\limits_{i=1}^{n}{{{x}_{i}}\ln {{y}_{i}}} \\ \sum\limits_{i=1}^{n}{\ln {{y}_{i}}} \\ \end{matrix} \right][/math]

Riešenie tejto sústavy rovníc je:

[math]\begin{align} & a=\frac{n\sum\limits_{i=1}^{n}{{{x}_{i}}\ln {{y}_{i}}}-\sum\limits_{i=1}^{n}{{{x}_{i}}}\sum\limits_{i=1}^{n}{\ln {{y}_{i}}}}{n\sum\limits_{i=1}^{n}{x_{i}^{2}}-{{\left( \sum\limits_{i=1}^{n}{{{x}_{i}}} \right)}^{2}}} \\ & B=\frac{\sum\limits_{i=1}^{n}{\ln {{y}_{i}}}-a\sum\limits_{i=1}^{n}{{{x}_{i}}}}{n} \\ \end{align}[/math]

kde [math]B=\ln b[/math].

Mocnimová aproximácia: [math]y=b{{x}^{a}}[/math]

Rovnicu [math]y=b{{x}^{a}}[/math] si upravíme tak, že ju zapíšeme v logaritmickom tvare:

[math]\begin{align} & \ln y=\ln \left( b{{x}^{a}} \right) \\ & \ln y=\ln b+\ln {{x}^{a}} \\ & \ln y=\ln b+a\ln x\,\,\left[ B=\ln b \right] \\ & \ln y=B+a\ln x \\ \end{align}[/math]

Dostávame sústavu 2 rovníc o 2 neznámych. V rovniciach je použitá substitúcia [math]B=\ln b[/math].

[math]\left[ \begin{matrix} \sum\limits_{i=1}^{n}{\ln {{x}_{i}}^{2}} & \sum\limits_{i=1}^{n}{\ln {{x}_{i}}} \\ \sum\limits_{i=1}^{n}{\ln {{x}_{i}}} & n \\ \end{matrix} \right]\cdot \left[ \begin{matrix} a \\ B \\ \end{matrix} \right]=\left[ \begin{matrix} \sum\limits_{i=1}^{n}{\ln {{x}_{i}}\ln {{y}_{i}}} \\ \sum\limits_{i=1}^{n}{\ln {{y}_{i}}} \\ \end{matrix} \right][/math]

Riešením je:

[math]\begin{align} & a=\frac{n\sum\limits_{i=1}^{n}{\ln {{x}_{i}}\ln {{y}_{i}}}-\sum\limits_{i=1}^{n}{\ln {{x}_{i}}}\sum\limits_{i=1}^{n}{\ln {{y}_{i}}}}{n\sum\limits_{i=1}^{n}{{{\left( \ln {{x}_{i}} \right)}^{2}}}-{{\left( \sum\limits_{i=1}^{n}{\ln {{x}_{i}}} \right)}^{2}}} \\ & B=\frac{\sum\limits_{i=1}^{n}{\ln {{y}_{i}}}-a\sum\limits_{i=1}^{n}{\ln {{x}_{i}}}}{n} \\ \end{align}[/math]

kde [math]B=\ln b[/math].

Analýza programátorského riešenia

Návrh dátových štruktúr

V úlohe budeme pracovať s bodmi v rovine, ktoré sú charakterizované x a y súradnicou. Preto si vytvoríme štruktúru Bod, ktorá bude takýto bod reprezentovať:

struct Bod{
   double x,y;
};

V zadaní je úlohou vypočítať rovnicu najlepšie aproximujúcej krivky (lineárnej, logaritmickej, exponenciálnej alebo mocninovej). Tieto krivky sú charakterizované (okrem svojej matematickej funkcie, ktorú vyjadríme neskôr) koeficientami a a b. Na to, aby sme rozlíšili o akú krivku ide, vytvorme si vymenovaný zoznam TYP_KRIVKY, ktorý bude ma 4 symbolické hodnoty:

enum TYP_KRIVKY{
   LIN, LOG, EXP, MOC
};

Význam týchto skratiek si zadefinujme nasledovne:

  • LIN - lineárna krivka (polynóm 1 stupňa) s rovnicou [math]y=ax+b[/math]
  • LOG - logaritmická krivka s rovnicou [math]y=a \ln{x} + b[/math]
  • EXP - exponenciálna krivka s rovnicou [math]y=b e^{ax}[/math]
  • MOC - mocninová krivka s rovnicou [math]y=bx^a[/math]

Teraz si môžeme zadefinovať štruktúru Krivka ako trojicu: typ (TYP_KRIVKY), koeficient a a koeficient b.

struct Krivka{
   TYP_KRIVKY typ;
   double a,b;
};

Riešenie v jazyku C

Vytvoríme funkciu Krivka minSq(Bod *body, int n) kde v poli body sú vstupné údaje, n je počet bodov. Funkcia bude vracať dátovú struktúru typu Krivka, kde budú informácie o tvare aproximujúcej krivky.

 1 Krivka minSq_lin(Bod *body, int n)
 2 {  // y=a+bx
 3    double sx=0,sxx=0,sy=0,sxy=0;
 4     for(int i=0;i<n;i++)
 5     {    sx+=body[i].x;
 6          sy+=body[i].y;
 7          sxy+=body[i].x * body[i].y;
 8          sxx+=body[i].x * body[i].x;
 9     }
10     Krivka k;
11     k.a=(n*sxy-sy*sx)/(n*sxx-sx*sx);
12     k.b=(sy-k.a*sx)/n;
13     k.typ=LIN;
14     return k;
15 }

Podobným spôsobom môžeme vytvoriť ďalšie 3 funkcie minSq_log, minSq_exp minSq_moc. Ukážeme si ale všeobecnejší tvar funkcie minSq, kde pridáme ďalší parameter, ktorý bude hovoriť aký typ aproximácie sa má vypočítať.

 1 Krivka minSq(Bod *body, int n, TYP_KRIVKY typ_proximacie)
 2 {  
 3    double sx=0,sxx=0,sy=0,sxy=0,slxly;
 4    double slx=0,slxx=0,sly=0,slxy=0,sxly=0;
 5     for(int i=0;i<n;i++)
 6     {    sx+=body[i].x;
 7          slx+=log(body[i].x);
 8          sy+=body[i].y;
 9          sly+=log(body[i].y);
10          sxy+=body[i].x * body[i].y;
11          slxy+=log(body[i].x) * body[i].y;
12          sxly+=body[i].x * log(body[i].y);
13          slxly+=log(body[i].x) *log(body[i].y);
14          sxx+=body[i].x * body[i].x;
15          slxx+=log(body[i].x) * log(body[i].x);
16     }
17     Krivka k;
18     switch(typ_proximacie)
19     {
20       case LIN:
21           k.a=(n*sxy-sy*sx)/(n*sxx-sx*sx);
22           k.b=(sy-k.a*sx)/n;
23           break;
24       case LOG:
25           k.a=(n*slxy-sy*slx)/(n*slxx-slx*slx);
26           k.b=(sy-k.a*slx)/n;
27           break;
28       case EXP:
29           k.a=(n*sxly-sly*sx)/(n*sxx-sx*sx);
30           k.b=(sly-k.a*sx)/n;
31           k.b=exp(k.b);
32           break;
33       case MOC:
34           k.a=(n*slxly-sly*slx)/(n*slxx-slx*slx);
35           k.b=(sly-k.a*slx)/n;
36           k.b=exp(k.b);
37           break;
38     }
39     k.typ=typ_proximacie;
40     return k;
41 }

Použitie v programe

Pre názornosť nebudeme v programe dáta načítavať z klávesnice, ale vstupné údaje budú zadané priamo v zdrojovom kóde. Modifikácia zdrojového kódu s načítaním vstupných hodnôt je jednoduchá a určite ju každý zvládne sám.

V zdrojovom kóde pribudla ešte funkcia vypisRovnicuKrivky, ktorá vypíše rovnicu krivky v čitateľnom tvare pre človeka. Vo výpise sú vynechané všetky funkcie, ktoré sú boli v tomto príklade vysvetlené a ich zdrojový kód bol napísaný skôr.

 1 #include <iostream.h>
 2 #include <conio.h>
 3 #include <math.h>
 4 
 5 struct Bod {double x,y;};
 6 enum TYP_KRIVKY{ LIN, LOG, EXP, MOC};
 7 struct Krivka{
 8    TYP_KRIVKY typ;
 9    double a,b;
10 };
11 
12 //Krivka minSq(Bod *body, int n, TYP_KRIVKY typ_proximacie)
13 
14 void vypisRovnicuKrivky(Krivka k)
15 {
16    cout<<"y=";   
17    switch(k.typ)
18    {
19       case LIN:
20            cout<<k.a<<" x";
21            if(k.b>0) cout<<"+";
22            cout<<k.b;
23           break;
24       case LOG:
25            cout<<k.a<<" ln x";
26            if(k.b>0) cout<<"+";
27            cout<<k.b;
28           break;
29       case EXP:
30            cout<<k.b<<" e^(";
31            cout<<k.a<<" x)";
32           break;
33       case MOC:
34            cout<<k.b<<" x^(";
35            cout<<k.a<<" )";
36           break;
37     }     
38 }
39 int main()
40 {  const int n=8;
41    Bod body[n]={
42       {1.50,0.58},
43       {2.00,8.61},
44       {3.00,60.91},
45       {4.00,201.40},
46       {5.00,476.59},
47       {6.00,892.57},
48       {7.00,1304.37},
49       {8.00,1115.04}		
50    }; 
51   Krivka Q;
52   Q=minSq(body,n,EXP); // LOG, LIN, EXP
53   vypisRovnicuKrivky(Q);
54   getch();
55 }

Overenie vypočítaných hodnôt

Na nasledujúcich obrázkoch sú vstupné body z príkladu znázornené v tabuľkovom procesore a je vypočítaná najlepšia aproximujúca krivka (MS Excel ju nazýva spojnica trendu).

Zadanie 2

Vytvorte funkciu, ktorá pre zadaná vstupy vypočíta všetky typy aproximácií (lineárnu, logaritmickú, exponenciálnu a mocninovú) a určí ktorá aproximácia je optimálna pre zadané vstupné body.

Analýza

Pri vyhodnocovaní optimálnej aproximácie musíme vypočítať všetky možné aproximačné krivky a zvoliť kritérium voľby najlepšej aproximácie. Kritériom posudzovania môže byť hľadanie minimálych hodnôt:

  • vzdialenosť uzlových (vstupných) bodov od aproximujúcej krivky
  • štvorce vzdialeností uzlových (vstupných) bodov od aproximujúcej krivky
  • plocha medziproximujúcou krivkou a interpolačným polynómom, ktorých prechádza cez uzlové body

Riešenie v jazyku C

Kriteriom na vyhodnotenie najlepšej aproximovanej krivky bude hľadanie minimálych hodnôt a vzdialeností uzlových (vstupných) bodov od aproximujúcej krivky.

Definovanie štruktúr potrebných pre spracovanie súradníc:

 1 struct Bod{
 2    double x,y;
 3 };
 4 
 5 enum TYP_KRIVKY{
 6    LIN, LOG, EXP, MOC
 7    };
 8 
 9 struct Krivka{
10    TYP_KRIVKY typ;
11    double a,b;
12 };

Táto štruktúra pre zjednodušenie poslúži na uloženie všetkých typov aproximácie do jednej premennej typu Aprox:

1 struct Aprox{
2     Krivka Lin, Log, Exp, Moc;
3 };

Funkcia aproximacie vypočíta všetky 4 typy aproximácií. Návratová hodnota funkcie je typu Aprox:

 1 Aprox aproximacie(Bod *body, int n, Aprox B)   //vypocita vsetky
 2                                                  //aproximacie
 3 {
 4    double sx=0,sxx=0,sy=0,sxy=0,slxly;
 5    double slx=0,slxx=0,sly=0,slxy=0,sxly=0;
 6 	for(int i=0;i<n;i++)
 7 	{    sx+=body[i].x;
 8 		 slx+=log(body[i].x);
 9 		 sy+=body[i].y;
10 		 sly+=log(body[i].y);
11 		 sxy+=body[i].x * body[i].y;
12 		 slxy+=log(body[i].x) * body[i].y;
13 		 sxly+=body[i].x * log(body[i].y);
14 		 slxly+=log(body[i].x) *log(body[i].y);
15          sxx+=body[i].x * body[i].x;
16          slxx+=log(body[i].x) * log(body[i].x);
17     }
18 
19 		 B.Lin.a=(n*sxy-sy*sx)/(n*sxx-sx*sx);
20 		 B.Lin.b=(sy-B.Lin.a*sx)/n;
21 
22 		 B.Log.a=(n*slxy-sy*slx)/(n*slxx-slx*slx);
23 		 B.Log.b=(sy-B.Log.a*slx)/n;
24 
25 		 B.Exp.a=(n*sxly-sly*sx)/(n*sxx-sx*sx);
26 		 B.Exp.b=(sly-B.Exp.a*sx)/n;
27 		 B.Exp.b=exp(B.Exp.b);
28 
29 		 B.Moc.a=(n*slxly-sly*slx)/(n*slxx-slx*slx);
30 		 B.Moc.b=(sly-B.Moc.a*slx)/n;
31 	     B.Moc.b=exp(B.Moc.b);
32 
33     return B;
34 }

Na vypísanie najvhodnejšieho typu rovnice je použitá funkcia vypisRovnicuKrivky:

 1 void vypisRovnicuKrivky(Krivka k)    //vypise rovnicu urciteho typu
 2 {
 3    cout<<"y = ";
 4    switch(k.typ)
 5    {
 6 	  case LIN:
 7 		   cout<<k.a<<" x";
 8 		   if(k.b>0) cout<<"+";
 9 		   cout<<k.b<<"   -> LINEARNA"<<endl;
10 		  break;
11 	  case LOG:
12 		   cout<<k.a<<" ln x";
13 		   if(k.b>0) cout<<"+";
14 		   cout<<k.b<<"   -> LOGARITMICKA"<<endl;
15 		  break;
16 	  case EXP:
17 		   cout<<k.b<<" e^(";
18 		   cout<<k.a<<" x)"<<"   -> EXPONENCIALNA"<<endl;
19 		  break;
20 	  case MOC:
21 		   cout<<k.b<<" x^(";
22 		   cout<<k.a<<" )"<<"   -> MOCNINOVA"<<endl;
23 		  break;
24 	}
25 }

Pre zistenie optimálnej aproximovanej krivky pre dané vstupné body je použitý jednoduchý algoritmus daný vo funkcii zistiOptimum.

  • Postup riešenia algoritmu spočíva v zistení vzdialenosti y-onovej súradnice vstupného bodu od aproximovanej y-onovej súradnice, ktorú musíme najskôr vypočítať, a porovnaním ako sú tieto body vzdialené od jednotlivých bodov kriviky pre každý typ, zistíme optimum.
  • Dosadením x-ových súradníc vstupných bodov do príslušnej rovnice krivky sa vypočíta aproximovaná y-onová súradnica. Absolútnu hodnotu rozdielu vzdialenosti medzi y-onovými súradnicami(vypočítanými a vstupnými) sčíta a uloží túto sumu do premennej sumLin (pre každý ďalší jeden typ aproximácie do - sumLog, sumExp, sumMoc ). Toto sa opakuje pre každý jeden vstupný bod a pre každú rovnicu krivky.
  • Na záver sa sumy všetkých typov aproximácií porovnajú. Najmenšia suma predstavuje najmenšie odchýlky pre danú krivku v porovnaní so vstupnými bodmi. Funkcia vráti Optimálnu krivku aproximácie s najmenšou sumou:
 1 Krivka zistiOptimum(Bod *body, int n, Aprox B)
 2 {
 3 	double sumLin=0, sumLog=0, sumExp=0, sumMoc=0, 	y1=0,y2=0,y3=0,y4=0, Naj=0;
 4 	Krivka Optim;
 5 	for(int i=0; i < n; i++)
 6 	{
 7 		y1=((B.Lin.a*body[i].x)+B.Lin.b)-body[i].y;  
 8 		if(y1<0) y1=-y1;
 9 		  sumLin+=y1;
10 
11 		y2=((B.Log.a*(log(body[i].x)))+B.Log.b)-body[i].y;  			
12                 if(y2<0) y2=-y2;                           
13 		  sumLog+=y2;
14 
15 		y3=(B.Exp.b*(exp(B.Exp.a*body[i].x)));     
16 		y3=log(y3);
17 		if(y3<0) y3=-y3;
18 		 sumExp+=y3;
19 
20 		y4=(B.Moc.b*(pow(body[i].x, B.Moc.a)))-body[i].y;
21 		if(y4<0) y4=-y4;                           
22 		  sumMoc+=y4;
23 
24 	}
25 	Naj=min(sumLin,sumLog);
26 	Naj=min(Naj,sumExp);
27 	Naj=min(Naj,sumMoc);
28 	
29 	if(Naj==sumLin) {B.Lin.typ=LIN; return B.Lin;}
30 	if(Naj==sumLog) {B.Log.typ=LOG; return B.Log;}
31 	if(Naj==sumExp) {B.Exp.typ=EXP; return B.Exp;}
32 
33 	B.Moc.typ=MOC;
34 	return B.Moc;
35 }

Modifikácia zdrojového kódu

Pre názornosť nebudeme v programe dáta načítavať z klávesnice, ale vstupné údaje budú zadané priamo v zdrojovom kóde

 1 #include <iostream.h>
 2 #include <conio.h>
 3 #include <math.h>
 4 #include <cmath>
 5 
 6 struct Bod{
 7    double x,y;
 8 };
 9 enum TYP_KRIVKY{
10    LIN, LOG, EXP, MOC
11    };
12 struct Krivka{
13    TYP_KRIVKY typ;
14    double a,b;
15 };
16 struct Aprox{
17 	Krivka Lin, Log, Exp, Moc;
18 };
19 
20 //Aprox aproximacie(Bod *body, int n, Aprox B)
21 //void vypisRovnicuKrivky(Krivka k)
22 //Krivka zistiOptimum(Bod *body, int n, Aprox B)
23 
24 int main(int argc, char* argv[])
25 {
26 	const int n=8;
27 	Bod body[n]={
28 	  {1.50,0.58},
29 	  {2.00,8.61},
30 	  {3.00,60.91},
31 	  {4.00,201.40},
32 	  {5.00,476.59},
33 	  {6.00,892.57},
34 	  {7.00,1304.37},
35 	  {8.00,2715.04}
36 	};
37 	Aprox B;
38 	cout<<"Najvhodnejsia: ";
39 	vypisRovnicuKrivky(zistiOptimum(body, n, B));
40 
41     getch();
42     return 0;
43 }