Numerické derivovanie (riešené príklady)

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

Vypočítajte hodnoty prvej derivácie 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:

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. 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:

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]

for(int i=0;i<4;i++)
 diferecie[i]=diferecie[i+1]-diferecie[i];

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]	[math]f(x_{0+2h})[/math]

Do vzorca pre deriváciu započítame 1. a 2. člen tohto poľa:

  double D;
  D=(diferecie[1]+diferecie[2])/2;

Druhá iterácia výpočtu diferencií: výpočet [math]{\Delta}^2 {y}_{i}[/math]

for(int i=0;i<3;i++)
 diferecie[i]=diferecie[i+1]-diferecie[i];

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]	[math]\Delta {y}_{2}[/math]	[math]f(x_{0+2h})[/math]

Tretia iterácia výpočtu diferencií: výpočet [math]{\Delta}^2 {y}_{i}[/math]

for(int i=0;i<2;i++)
 diferecie[i]=diferecie[i+1]-diferecie[i];

index	0	1	2	3	4
hodnota	[math]{\Delta}^3 {y}_{-1}[/math]	[math]{\Delta}^3 {y}_{1}[/math]	[math]{\Delta}^2 {y}_{1}[/math]	[math]\Delta {y}_{2}[/math]	[math]f(x_{0+2h})[/math]

Do vzorca pre deriváciu započítame 0. a 1. člen tohto poľa:

  D-=(diferecie[0]-diferecie[1])/12;

Dopočítame konečnú hodnotu derivácie v bode x₀:

  D=D/h;

Výsledný funkcia môže mať tvar:

 1 double derivaciaL(Bod *data, int n,double x0, double h=0.0001)
 2 {
 3   double diferecie[5],D=0;
 4   int i;
 5   for(i=-2;i<=2;i++)
 6     diferecie[i+2]=LagrangeInterpol(data,n,x0+i*h);
 7   for(i=0;i<4;i++)
 8     diferecie[i]=diferecie[i+1]-diferecie[i];
 9   D=(diferecie[1]+diferecie[2])/2;
10   for(i=0;i<3;i++)
11     diferecie[i]=diferecie[i+1]-diferecie[i];
12   for(i=0;i<2;i++)
13     diferecie[i]=diferecie[i+1]-diferecie[i];
14   D-=(diferecie[0]-diferecie[1])/12;
15   return D/h;
16 }

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

[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.

 1 #include <iostream.h>
 2 #include <conio.h>
 3 #include <math.h>
 4 struct Bod {double x,y;};
 5 
 6 //double LagrangeInterpol(Bod *data, int n, double x)
 7 //double NewtonPol(Bod *data, int n,double x)
 8 
 9 double derivacia(double (*polynom)(Bod*,int,double),Bod *data, int n,double x0, double h=0.01)
10 {
11   double diferecie[5],D=0;
12   int i;
13   for(i=-2;i<=2;i++)
14    { diferecie[i+2]=polynom(data,n,(x0+i*h));
15      //cout<<diferecie[i];
16    }
17   for(i=0;i<4;i++)
18     diferecie[i]=diferecie[i+1]-diferecie[i];
19   D=(diferecie[1]+diferecie[2])/2;
20   for(i=0;i<3;i++)
21     diferecie[i]=diferecie[i+1]-diferecie[i];
22   for(i=0;i<2;i++)
23     diferecie[i]=diferecie[i+1]-diferecie[i];
24   D-=(diferecie[0]+diferecie[1])/12;
25   return D/h;
26 }
27 int main()
28 {  const int n=8;
29    Bod body[n]={
30       {1.0, -1.71828},
31       {1.3, -0.8132},
32       {2.1, 11.28193},
33       {3.8, 163.8124},
34       {4.55, 333.9611},
35       {5.0, 476.5868},
36       {6.7, 1202.706},
37       {7.9 , 1197.726} }; 
38   double x,krok=0.2;
39  
40   for (x=body[0].x ; x<body[7].x ; x+=krok)
41   {
42       cout<<"Newton. polynom: f'("<<x<<")="<<derivacia(NewtonPol,body,n,x)<<endl;
43       cout<<"Lagran. polynom: f'("<<x<<")="<<derivacia(LagrangeInterpol,body,n,x)<<endl;
44       cout<<"\trozdiel:"<<derivacia(NewtonPol,body,n,x)-derivacia(LagrangeInterpol,body,n,x)<<endl;
45   }
46   getch();
47 }