Metóda najmenších štvorcov pre aproximáciu dát pomocou systému Mathematica

Základy

Dané je: Tabuľka funkcie  y= f(x): { (Subscript[x, 1]),(Subscript[y, 1])}, ...,{ (Subscript[x, n1]),(Subscript[y, n1])} .

Úloha: Nájdite algebrický polynóm daného  stupňa, 
ktorý najlepšie aproximuje uvedenú funkciu,
použitím metódy najmenších štvorcov (LSM)určeného stupňa.

Krátky úvod do teórie:
Metódu navrhol Carl Friedrich Gauss v roku 1794.  
Je známa aj pod názvom metóda  regresnej analýzy.

Hľadaný polynóm  má tvar:  (Subscript[p, k])(x)= (Subscript[c, 0]) +(Subscript[c, 1])x +(Subscript[c, 2])(x^2)+...+ (Subscript[c, k])(x^k),
pričom koeficienty  (Subscript[c, 0]) , (Subscript[c, 1]),..., (Subscript[c, k])  vyhovujú sústave:
S.c=d,
S=(({{(Subscript[s, 0]), (Subscript[s, 1]), (Subscript[s, 2]), ., ., (Subscript[s, k])}, {(Subscript[s, 1]), (Subscript[s, 2]), (Subscript[s, 3]), ., ., (Subscript[s, k+1])}, {(Subscript[s, 2]), (Subscript[s, 3]), (Subscript[s, 4]), ., ., (Subscript[s, k+2])}, {., ., ., ., ., .}, {., ., ., ., ., .}, {(Subscript[s, k]), (Subscript[s, k+1]), ., ., ., (Subscript[s, 2k])}})),   c=(({{(Subscript[c, 0])}, {(Subscript[c, 1])}, {.}, {.}, {.}, {(Subscript[c, k])}})), d=(({{(Subscript[d, 0])}, {(Subscript[d, 1])}, {.}, {.}, {.}, {(Subscript[d, k])}})),   
    (Subscript[s, m])=  (∑_(i=1)^n1)((Subscript[x, [[i]]])^m),  m=0,...,2k;        (Subscript[d, m])=(∑_(i=1)^n1)((Subscript[x, [[i]]])^m)*(Subscript[y, [[i]]]),  m=0,...,k.

Stredná kvadratická odchýlka je:    (Subscript[δ, k])=(Sqrt[(∑_(i=1)^n1)(((f((Subscript[x, i]))-(Subscript[p, k])((Subscript[x, i]))))^2)])  .

Aplikácie: Polynóm p(x) sa používa na aproximáciu  funkcie f(x) 
v okolí nejakého bodu  t  tak, že platí  (Subscript[p, k])(t)≈f(t).

Príklad 1.    Nech je    f(x)=x+sin((x^2)), 1≤x≤2.   
a) Vytvorte tabuľku hodnôt funkcie pre 11 rovnomerne rozložených bodov intervalu.
 b)  Nájdite polynóm prvého stupňa, ktorý najlepšie aproximuje funkciu vzhľadom na LSM.
        c) Vyhodnoďte chybu aproximácie.
               d) Vypočítajte približnú hodnotu funkcie  f(x)  pre  x=1.33 , pomocou nájdeného polynómu.

Riešenie:
Definujeme funkciu a nakreslíme jej graf na intervale [1,2]:

f[xx_]:=xx+Sin[(xx^2)]    
a=1.;b=2. ;
gf=Plot[f[xx],{xx,a,b}]

a)Vypočítame hodnoty členov postuponosti bodov :  (Subscript[x, i]),(Subscript[y, i])=f((Subscript[x, i])), i=1,2,...n1, pre zvolených  n1=11 rovnomerne rozložených bodov s krokom  h na intervale  [a,b]:

n1=11;  h=(((b-a))/(n1-1));Print["Step h= ",h]
x=Table[a+h*(i-1), {i,1,n1}]; Print["x=",x] ;
y=Table[f[(Subscript[x, [[i]]])], {i,1,n1}]; Print["y=",y] ;
xy=Table[{(Subscript[x, [[i]]]),(Subscript[y, [[i]]])},{i,1,n1}]
gy=ListPlot[xy,PlotStyle->PointSize[0.02]]

b)Vypočítame všetky súčty potrebné pre LSM prvého stupňa, príslušný systém vyriešime pomocou príkazu LinearSolve[ ] v systéme Mathematica  a extrahujeme hodnoty koeficientov polynómu - c0 a c1:

s0=n1; s1=(∑_(i=1)^n1)(Subscript[x, [[i]]])  ;s2=(∑_(i=1)^n1)((Subscript[x, [[i]]])^2);d0=(∑_(i=1)^n1)(Subscript[y, [[i]]]);d1=(∑_(i=1)^n1)(Subscript[x, [[i]]])*(Subscript[y, [[i]]]);
S=(({{s0, s1}, {s1, s2}}))
d=(({{d0}, {d1}}))
c=LinearSolve[S,d]
c0=(Subscript[(Subscript[c, [[1]]]), [[1]]])
c1=(Subscript[(Subscript[c, [[2]]]), [[1]]])

Definujeme polynóm, ktorý je najlepšou  aproximáciou danej funkcie pomocou metódy LSM prvého stupňa
 a nakreslíme jeho graf. 
Napokon použijeme príkaz Show[]  a zobrazíme oba grafy v jednom obrázku,  teda graf tabuľkovej funkcie   f(x)  a súčasne graf jej aproximácie  p1(x):

p1[xx_]:=c0+c1*xx
gp1=Plot[p1[xx],{xx,a,b},PlotStyle->Hue[.6]]
Show[gy,gp1]

c) Vypočítame strednú kvadratickú odchýlku ako chybu aproximácie:

(Subscript[δ, 1])=(Sqrt[(∑_(i=1)^n1)(((f[(Subscript[x, [[i]]])]-p1[(Subscript[x, [[i]]])]))^2)])

d) Takto môžeme vypočítať približnú hodnotu funkcie v ľubovoľnom bode   f(1.33). V mnohých konkrétnych prípadoch nepoznáme presný funkčný predpis a musíme sa spoľahnúť iba na tabuľku jej nameraných hodnôt:

p1[1.33]  
Abs[f[1.33]-p1[1.33] ]

Záver:   Vypočítaná chyba (Subscript[δ, 1]) je evidentne veľká a je zrejmá aj zo zobrazených grafov oboch funkcií. Je preto potrebné nájsť aproximáciu polynómom vyššieho stupňa.

Príklad 2.  Vrátime sa k funkcii z predchádzajúceho príkladu 1:
   a)  Nájdeme polynóm druhého stupňa, ktorý je najlepšou  aproximáciou pomocou LSM. 
b) Vyhodnotíme chybu aproximácie. 
c) Nájdeme približnú hodnotu  f(1.33).

Riešenie :

a) Vypočítame všetky súčty potrebné pre metódu LSM druhého stupňa, vyriešime príslušný systém pomocou príkazu LinearSolve a nájdeme koeficienty polynómu, ktoré označíme v0, v1 a v2:

s0=n1; s1=(∑_(i=1)^n1)(Subscript[x, [[i]]])  ;s2=(∑_(i=1)^n1)((Subscript[x, [[i]]])^2);s3=(∑_(i=1)^n1)((Subscript[x, [[i]]])^3);s4=(∑_(i=1)^n1)((Subscript[x, [[i]]])^4);
d0=(∑_(i=1)^n1)(Subscript[y, [[i]]]);d1=(∑_(i=1)^n1)(Subscript[x, [[i]]])*(Subscript[y, [[i]]]);d2=(∑_(i=1)^n1)((Subscript[x, [[i]]])^2)*(Subscript[y, [[i]]]);
S=(({{s0, s1, s2}, {s1, s2, s3}, {s2, s3, s4}}))
d=(({{d0}, {d1}, {d2}}))
c=LinearSolve[S,d]
v0=(Subscript[(Subscript[c, [[1]]]), [[1]]])
v1=(Subscript[(Subscript[c, [[2]]]), [[1]]])
v2=(Subscript[(Subscript[c, [[3]]]), [[1]]])

Definujeme polynóm druhého stupňa, ktorý je najlepšou aproximáciou funkcie vzhľadom na LSM,a nakreslíme jeho graf. Napokon použijeme príkaz Show[]  a nakreslíme súčasne oba grafy, graf tabuľkovej funkcie  f(x) a jej  aproximácie p2(x):

p2[xx_]:=v0+v1*xx+v2*(xx^2)
gp2=Plot[p2[xx],{xx,a,b},PlotStyle->Hue[.3]]
Show[gy,gp2]

b) Vypočítame strednú kvadratickú odchýlku - chybu aproximácie:

(Subscript[δ, 2])=(Sqrt[(∑_(i=1)^n1)(((f[(Subscript[x, [[i]]])]-p2[(Subscript[x, [[i]]])]))^2)])

c) Nájdeme približnú hodnotu funkcie  f(1.33) a absolútnu chybu :

p2[1.33]  
aer=Abs[f[1.33]-p2[1.33]]; Print[" The test absolute error now is=", aer]

Záver:   Získané chyby  (Subscript[δ, 2]) a aer asú evidentne ešte stále veľké, ale akceptovateľnejšie, ako v predchádzajúcej  aproximácii (porovnajte aj grafy oboch funkcií).  Nájdeme preto aproximáciu polynómom ešte vyššieho stupňa.

Príklad 3.  Opäť sa vrátime k funkcii z príkladu 1:
   a) Nájdeme aproximačný polynóm piateho stupňa .
b) Vyhodnotíme chybu  aproximácie.  
c) Nájdeme približnú hodnotu funkcie  f(1.33).

Riešenie:

a) Vypočítanme všetky potrebné súčty pre LSM piateho stupňa, vyriešime algebrický lineárny systém pomocou príkazu LinearSolve[] 
a extrahujeme hodnoty koeficientov polynómu, ktoré označíme  z0, z1, z2 , z3,z4 a z5:

s0=n1; s1=(∑_(i=1)^n1)(Subscript[x, [[i]]])  ;s2=(∑_(i=1)^n1)((Subscript[x, [[i]]])^2);s3=(∑_(i=1)^n1)((Subscript[x, [[i]]])^3);s4=(∑_(i=1)^n1)((Subscript[x, [[i]]])^4);
s5=(∑_(i=1)^n1)((Subscript[x, [[i]]])^5);s6=(∑_(i=1)^n1)((Subscript[x, [[i]]])^6);s7=(∑_(i=1)^n1)((Subscript[x, [[i]]])^7);
s8=(∑_(i=1)^n1)((Subscript[x, [[i]]])^8);s9=(∑_(i=1)^n1)((Subscript[x, [[i]]])^9);s10=(∑_(i=1)^n1)((Subscript[x, [[i]]])^10);
d0=(∑_(i=1)^n1)(Subscript[y, [[i]]]);d1=(∑_(i=1)^n1)(Subscript[x, [[i]]])*(Subscript[y, [[i]]]);d2=(∑_(i=1)^n1)((Subscript[x, [[i]]])^2)*(Subscript[y, [[i]]]);
d3=(∑_(i=1)^n1)((Subscript[x, [[i]]])^3)*(Subscript[y, [[i]]]);d4=(∑_(i=1)^n1)((Subscript[x, [[i]]])^4)*(Subscript[y, [[i]]]);d5=(∑_(i=1)^n1)((Subscript[x, [[i]]])^5)*(Subscript[y, [[i]]]);
S=(({{s0, s1, s2, s3, s4, s5}, {s1, s2, s3, s4, s5, s6}, {s2, s3, s4, s5, s6, s7}, {s3, s4, s5, s6, s7, s8}, {s4, s5, s6, s7, s8, s9}, {s5, s6, s7, s8, s9, s10}}));d=(({{d0}, {d1}, {d2}, {d3}, {d4}, {d5}}));
c=LinearSolve[S,d]
z0=(Subscript[(Subscript[c, [[1]]]), [[1]]]);z1=(Subscript[(Subscript[c, [[2]]]), [[1]]]);z2=(Subscript[(Subscript[c, [[3]]]), [[1]]]);z3=(Subscript[(Subscript[c, [[4]]]), [[1]]]);
z4=(Subscript[(Subscript[c, [[5]]]), [[1]]]);z5=(Subscript[(Subscript[c, [[6]]]), [[1]]]);

Definujeme najlepšiu aproximáciu polynómom piateho stupňa a nakreslíme graf. 
Použijeme príkaz Show[]  a zobrazíme súčasne graf funkcie   f(x) a jej  aproximácie p5(x):

p5[xx_]:=z0+z1*xx+z2*(xx^2)+z3*(xx^3)+z4*(xx^4)+z5*(xx^5)
gp5=Plot[p5[xx],{xx,a,b},PlotStyle->Hue[.8]]
Show[gy,gp5]

b) Vypočítame strednú kvadratickú odchýlku aproximácie:

(Subscript[δ, 5])=(Sqrt[(∑_(i=1)^n1)(((f[(Subscript[x, [[i]]])]-p5[(Subscript[x, [[i]]])]))^2)])

c) Približná hodnota funkcie   f(1.33) je:

p5[1.33]  
aer5=Abs[f[1.33]-p5[1.33]]; Print[" The test absolute error now is=", aer5]

Záver:   Získaná chyba (Subscript[δ, 5]) a absolútna chyba majú pomerne vyhovujúce hodnoty.

Príklad 4.  Zobrazte riešenia predchádzajúcich troch príkladov v jednom grafe.

Riešenie:

Show[gy,gp1,gp2,gp5]