Lab. 3 Różniczkowanie numeryczne

Definicja pochodnej

$$f^{\prime}(x)=\lim _{h \rightarrow 0} \frac{f(x+h)-f(x)}{h}$$

Z szeregu Tylora

$$f(x+h)=f(x)+h f^{\prime}(x)+\frac{h^{2}}{2 !} f^{\prime \prime}(x)+\frac{h^{3}}{3 !} f^{\prime \prime \prime}\left(\xi_{1}\right)$$

Stąd przybliżenie

$$f^{\prime}(x)=\frac{f(x+h)-f(x)}{h} + O(h)$$

Dwupunktowe różnice

• dwupunktowe różnice zwykłe

$$f^{\prime}\left(x_{i}\right)=\frac{f\left(x_{i+1}\right)-f\left(x_{i}\right)}{h} + O(h)$$

• dwupunktowe różnice wsteczne

$$f^{\prime}\left(x_{i}\right)=\frac{f\left(x_{i}\right)-f\left(x_{i-1}\right)}{h} + O(h)$$

• dwupunktowe różnice centralne

$$f^{\prime}\left(x_{i}\right)=\frac{f\left(x_{i+1}\right)-f\left(x_{i-1}\right)}{2 h} + O(h^2)$$

Trzypunktowe różnice

• trzypunktowe różnice zwykłe

$$f^{\prime}\left(x_{i}\right)=\frac{-3 f\left(x_{i}\right)+4 f\left(x_{i+1}\right)-f\left(x_{i+2}\right)}{2 h} + O(h^2)$$

• trzypunktowe różnice wsteczne

$$f^{\prime}\left(x_{i}\right)=\frac{f\left(x_{i-2}\right)-4 f\left(x_{i-1}\right)+ 3 f\left(x_{i}\right)}{2 h} + O(h^2) $$

Przykład

Pochodna funkcji $f(x) = e^x$

Wartość dokładna $f'(x) = e^x$

exp_diff.c

#include <stdio.h>
#include <math.h>
 
float exp_diff_for(float x, float h)
{
    return (exp(x+h) - exp(x))/h ;
}
 
double exp_diff_back(double x, double h)
{
    return (exp(x-h) - exp(x))/h ;
}
 
double exp_diff_central(double x, double h)
{
    return (exp(x+h) - exp(x-h))/(2.0 * h);
}
 
int main()
{
  int n;
  float x = 1.0 , h, fx, err, ex;
 
  ex = exp(x);
 
  printf("# exact = %f\n", ex);
  printf("#   h         diff            error\n");
  for(h = 1e-1; h > 1e-6; h = h / 2.0)
  {
      fx = exp_diff_for(x, h);
      err = fabs((ex - fx) / ex);
      printf("%10.7lf %15.10lf %15.10lf\n", h, fx, err);
  }
 
  return 0;
}

Błąd zaokrąglenia i obcięcia

Druga pochodna

• różnice trzypunktowe zwykłe

$$f^{\prime \prime}\left(x_{i}\right)=\frac{f\left(x_{i}\right)-2 f\left(x_{i+1}\right)+f\left(x_{i+2}\right)}{h^{2}} + O(h)$$

• różnice trzypunktowe wsteczne

$$f^{\prime \prime}\left(x_{i}\right)=\frac{f\left(x_{i-2}\right)-2 f\left(x_{i-1}\right)+f\left(x_{i}\right)}{h^{2}} + O(h)$$

• różnice trzypunktowe centralne

$$f^{\prime \prime}\left(x_{i}\right)=\frac{f\left(x_{i-1}\right)-2 f\left(x_{i}\right)+ f\left(x_{i+1}\right)}{h^{2}} + O(h^2)$$

Zadanie:

Plik rocket.txt zawiera pomiary wysokości rakiety (w km) od momentu wystrzelenia t=0 s. do t=200 s. Plik zawiera pomiary zapisane w następującym formacie:

0 0
20 9.170
40 23.835

Każda linia zwiera pomiar czasu $t_i$ oraz wysokość rakiety.

Napisz program, który wczyta dane z pliku a następnie wyznaczy numerycznie wartości prędkości (pierwsza pochodna) oraz przyśpieszenia (druga pochodna) w każdym punkcie czasu $t_i = 0, 20, 40, \ldots$

Do wyznaczenia pochodnej użyj wzorów o dokładności co najmniej $O(h^2)$ (np. trzypunktowe różnice). Pochodną dla początkowych wartości $t_i=0$ wyznacz z różnic zwykłych Wartości pochodnej dla ostatniego elementu ($t_i=200$ s.) wyznacz z różnic wstecznych. Pozostałe wartości wyznacz z różnic centralnych.

Operacje na plikach w C

read.c

#include <stdio.h>
#define MAX 100
 
int main()
{
  int i, k;
  float tab[MAX];
  FILE *plik;
 
  plik = fopen("text.txt", "r");
  if ( plik == NULL) 
  {
     printf("Blad otworzenia pliku.");
     return 1;
  }
 
  i = 0;
  while(!feof(plik) && i < MAX)
     if(fscanf(plik, "%f", tab+i) > 0) i++;
 
  for(k=0; k<i; k++) 
     printf("%i -> %f\n", k, tab[k]);
 
  fclose(plik);
 
  return 0;
}

Dodatkowe ćwiczenia

Ćw. 1

Policzyć pochodną funkcji $f(X) = x \sin{x} $ ⁡w punkcie $x= \pi$ korzystając z metody dwupunktowej „w przód” i „w tył” oraz różnicy centralnej. Zbadaj jak zmienia się błąd względny i bezwzględny wyznaczenia pochodnej dla różnych wartości $h = 1, 0.1, 0.01, 0.001, 0.0001, 0.00001$.

Ćw. 2

Policzyć drugą pochodną funkcji $f(X) = x \sin{x} $ ⁡w punkcie $x= \pi$ korzystając z metody trzypunktowej centralnej. Zbadaj jak zmienia się błąd względny i bezwzględny wyznaczenia pochodnej dla różnych wartości $h = 1, 0.1, 0.01, 0.001, 0.0001, 0.00001$.

Ćw. 3

Plik heart_rate.csv zawiera przykład elektrokardiogramu. Sygnał jest rejestrowany z częstotliwością 100Hz. Napisz program, który będzie wykrywał uderzenia serca na podstawie pochodnej napięcia mierzonego na elektrodzie a następnie wyznaczy średnią wartość uderzeń na minutę.