Zobrazenie zdroja stránky Numerické derivovanie (riešené príklady)

{{Draft}}
{{Skripta programovanie (zbierka úloh)}}

==Zadanie==
Vypočítajte hodnoty prvej derivácie [[Algoritmy numerickej interpolácie (riešené príklady)|interpolačného polynómu]] v Newtonovom a Lagrangeovom tvare. Porovnajte tieto hodnoty.
Pre výpočet derivácie použite metódu rozdielových diferencií.

==Analýza==
Pre výpočet prvej derivácie pomocou rozdielových diferencií je vzťah:
:<math>{f}'\left( {{x}_{0}} \right)=\frac{1}{h}\left( \frac{\Delta {{y}_{-1}}+\Delta {{y}_{1}}}{2}-\frac{1}{6}\frac{{{\Delta }^{3}}{{y}_{-1}}+{{\Delta }^{3}}{{y}_{1}}}{2} \right)</math>

kde <math>\Delta {{y}_{-1}}=f(x_0-h)-f(x_0)</math>, 

<math>\Delta {{y}_{1}}=f(x_0)-f(x_0+h)</math>, 

<math>{{\Delta }^{3}}{{y}_{-1}}={{\Delta }^{2}}{{y}_{1}}-{{\Delta }^{2}}{{y}_{0}}</math>, 

<math>{{\Delta }^{3}}{{y}_{1}}={{\Delta }^{2}}{{y}_{0}}-{{\Delta }^{2}}{{y}_{1}}</math>.

Celá schéma výpočtu diferencií je na nasledujúcom obrázku:

[[Súbor:prog deriv diff.png|center|framed|Princíp výpočtu diferencií]]

Interpolačný polynóm v Lagrangeovom a Newtonovom tvare sú opísané a implementované v kapitole [[Algoritmy numerickej interpolácie (riešené príklady)|Algoritmy numerickej interpolácie]]. V tomto príklade použijeme už vytvorené funkcie ''LagrangeInterpol'' a ''NewtonPol''.

==Riešenie v jazyku C==
Vytvoríme si funkciu ''derivacia'', ktorej parametrami budú vstupné body pre interpoláciu a bod v ktorom chceme deriváciu počítať. Vo funkcii vypočítame hodnotu interpolačného polynómu pre potrebné body na osi x a následne vypočítame prvú deriváciu.
===Výpočet rozdielových diferencií===
Pre daný bod <math>x_0</math> potrebujeme vedieť funkčnú hodnotu v tomto bode a následne funkčné hodnoty o vzdialenosť h a 2h od tohto bodu <math>x_0</math> do ľavej aj pravej strany. Teda spolu 5 hodnôt. Vytvoríme si teda pomocné pole reálnych čísel o veľkosti 5. Do tohoto poľa vložíme funkčné hodnoty interpolačného polynómu:
{| class=prettytable
|+ double diferecie[5];
|-
!index
|0
|1
|2
|3
|4
|-
!hodnota
|<math>f(x_{0-2h})</math>
|<math>f(x_{0-h})</math>
|<math>f(x_{0})</math>
|<math>f(x_{0+h})</math>
|<math>f(x_{0+2h})</math>
|}
Poznámka: funkčné hodnoty f(x) dostaneme z interpolačného polynómu. Parameter h volíme malý (blízky 0).

'''Prvá iterácia''' výpočtu diferencií: výpočet <math>\Delta {y}_{i}</math>
<source lang="c">
for(int i=0;i<4;i++)
 diferecie[i]=diferecie[i+1]-diferecie[i];
</source>
{| class=prettytable
|-
!index
|0
|1
|2
|3
|4
|-
!hodnota
|<math>\Delta {y}_{-2}</math>
|<math>\Delta {y}_{-1}</math>
|<math>\Delta {y}_{1}</math>
|<math>\Delta {y}_{2}</math>
|style="color:gray"|<math>f(x_{0+2h})</math> 
|}
Do vzorca pre deriváciu započítame 1. a 2. člen tohto poľa:
<source lang="c">
  double D;
  D=(diferecie[1]+diferecie[2])/2;
</source>

'''Druhá iterácia''' výpočtu diferencií: výpočet <math>{\Delta}^2 {y}_{i}</math>
<source lang="c">
for(int i=0;i<3;i++)
 diferecie[i]=diferecie[i+1]-diferecie[i];
</source>
{| class=prettytable
|-
!index
|0
|1
|2
|3
|4
|-
!hodnota
|<math>{\Delta}^2 {y}_{-1}</math>
|<math>{\Delta}^2 {y}_{0}</math>
|<math>{\Delta}^2 {y}_{1}</math>
|style="color:gray"|<math>\Delta {y}_{2}</math>
|style="color:gray"|<math>f(x_{0+2h})</math> 
|}
'''Tretia iterácia''' výpočtu diferencií: výpočet <math>{\Delta}^2 {y}_{i}</math>
<source lang="c">
for(int i=0;i<2;i++)
 diferecie[i]=diferecie[i+1]-diferecie[i];
</source>
{| class=prettytable
|-
!index
|0
|1
|2
|3
|4
|-
!hodnota
|<math>{\Delta}^3 {y}_{-1}</math>
|<math>{\Delta}^3 {y}_{1}</math>
|style="color:gray"|<math>{\Delta}^2 {y}_{1}</math>
|style="color:gray"|<math>\Delta {y}_{2}</math>
|style="color:gray"|<math>f(x_{0+2h})</math> 
|}
Do vzorca pre deriváciu započítame 0. a 1. člen tohto poľa:
<source lang="c">
  D-=(diferecie[0]-diferecie[1])/12;
</source>
Dopočítame konečnú hodnotu derivácie v bode x<sub>0</sub>:
<source lang="c">
  D=D/h;
</source>

Výsledný funkcia môže mať tvar:
<source lang="c" line>
double derivaciaL(Bod *data, int n,double x0, double h=0.0001)
{
  double diferecie[5],D=0;
  int i;
  for(i=-2;i<=2;i++)
    diferecie[i+2]=LagrangeInterpol(data,n,x0+i*h);
  for(i=0;i<4;i++)
    diferecie[i]=diferecie[i+1]-diferecie[i];
  D=(diferecie[1]+diferecie[2])/2;
  for(i=0;i<3;i++)
    diferecie[i]=diferecie[i+1]-diferecie[i];
  for(i=0;i<2;i++)
    diferecie[i]=diferecie[i+1]-diferecie[i];
  D-=(diferecie[0]-diferecie[1])/12;
  return D/h;
}
</source>
Funkciu derivaciaN vytvoríme rovnako, ale na riadku č 6. bude volanie funkcie pre výpočet interpolačného polynómu v Newtonovom tvare.

==Vstupné údaje==
Pre potreby porovnania budú vstupné údaje rovnaké ako pri implementácii interpolačných polynómov:
Pre účel zadania použime nasledujúce vstupné hodnoty:
n = 8
{| class=datatable
|-
!<math>x_i</math>
!<math>y_i</math>
|-
|1.0
| -1.71828
|-
|1.3
| -0.8132
|-
|2.1
|11.28193
|-
|3.8
|163.8124
|-
|4.55
|333.9611
|-
|5.0
|476.5868
|-
|6.7
|1202.706
|-
|7.9
|1197.726
|}
Poznámka: vstupné body patria funkcii <math>4x^3-e^x</math>

==Celkové riešenie v jazyku C==
Budeme počítať hodnoty derivácie interpolačného polynómu na intervale určenom krajnými vstupnými bodmi s krokom 0.2.
<source lang="c" line>
#include <iostream.h>
#include <conio.h>
#include <math.h>
struct Bod {double x,y;};

//double LagrangeInterpol(Bod *data, int n, double x)
//double NewtonPol(Bod *data, int n,double x)

double derivacia(double (*polynom)(Bod*,int,double),Bod *data, int n,double x0, double h=0.01)
{
  double diferecie[5],D=0;
  int i;
  for(i=-2;i<=2;i++)
   { diferecie[i+2]=polynom(data,n,(x0+i*h));
     //cout<<diferecie[i];
   }
  for(i=0;i<4;i++)
    diferecie[i]=diferecie[i+1]-diferecie[i];
  D=(diferecie[1]+diferecie[2])/2;
  for(i=0;i<3;i++)
    diferecie[i]=diferecie[i+1]-diferecie[i];
  for(i=0;i<2;i++)
    diferecie[i]=diferecie[i+1]-diferecie[i];
  D-=(diferecie[0]+diferecie[1])/12;
  return D/h;
}
int main()
{  const int n=8;
   Bod body[n]={
      {1.0, -1.71828},
      {1.3, -0.8132},
      {2.1, 11.28193},
      {3.8, 163.8124},
      {4.55, 333.9611},
      {5.0, 476.5868},
      {6.7, 1202.706},
      {7.9 , 1197.726} }; 
  double x,krok=0.2;
 
  for (x=body[0].x ; x<body[7].x ; x+=krok)
  {
      cout<<"Newton. polynom: f'("<<x<<")="<<derivacia(NewtonPol,body,n,x)<<endl;
      cout<<"Lagran. polynom: f'("<<x<<")="<<derivacia(LagrangeInterpol,body,n,x)<<endl;
      cout<<"\trozdiel:"<<derivacia(NewtonPol,body,n,x)-derivacia(LagrangeInterpol,body,n,x)<<endl;
  }
  getch();
}
</source>