Lab. 2 Błędy numeryczne

Rodzaje błędów numerycznych

• błędy zaokrąglenia Round-off_error wynikające ze skończonej precyzji obliczeń
• błędy obcięcia Truncation_error wynikające z przybliżenia stosowanego w obliczeniach matematycznych

Utrata cyfr znaczących Loss of significance - znaczny wzrost bledu względnego w wyniku operacji na liczbach w skończonej precyzji

• błędy mogą się akumulować w trakcie obliczeń
• źle uwarunkowane zadania → nieograniczony wzrost błędu

Błąd względny i bezwzględny

Błąd bezwzględny $$\Delta x=\left|x-x_{0}\right|$$

Błąd względny $$\delta=\frac{\Delta x}{x}=\frac{\left|x-x_{0}\right|}{x}$$

możliwe do wyznaczenia tylko gdy znamy wartość dokładną $x$

Obliczanie wartości funkcji

Funkcja $e^{-x}$

$$e^{-x}= 1 - x + \frac{x^{2}}{2} - \frac{x^{3}}{6} + \cdots = \sum_{n=0}^{\infty}(-1)^{n} \frac{x^{n}}{n !}$$

Przykład: $e^{-x}$ metodą wprost

exp_direct.c

#include <stdio.h>
#include <math.h>

#define TRUNCATION 1e-10

double factorial(int n)
{
  int loop;
  double fac;
  for (loop = 1, fac = 1.0; loop <= n; loop++)
  {
    fac *= loop;
  }
  return fac;
}

int main()
{
  int n;
  double x, term, sum;

  printf("       x       exp(-x)       series       terms\n");

  for (x = 0.0; x < 100.0; x += 10.0)
  {
    sum = 0.0;

    n = 0;
    term = 1;
    while (fabs(term) > TRUNCATION)
    {
      term = pow(-1, n) * (pow(x, n) / factorial(n));
      sum += term;
      n++;
    }

    printf("%12f %12g %12g %5d\n", x, exp(-x), sum, n - 1);
  }
  return 0;
}

Przykład: $e^{-x}$ metodą rekurencyjną

$$e^{-x}=\sum_{n=0}^{\infty}(-1)^{n} \frac{x^{n}}{n !}=\sum_{n=0}^{\infty} s_{n} $$ gdzie $$ s_{n}=-s_{n-1} \frac{x}{n} $$ $$\quad s_0 = 1$$

exp_recursive.c

#include <stdio.h>
#include <math.h>

#define TRUNCATION 1e-10

int main()
{
  int loop, n;
  double x, term, sum;

  printf("       x       exp(-x)       series       terms\n");

  for (loop = 0; loop <= 100; loop += 10)
  {
    x = (double)loop;

    sum = 1.0;
    term = 1;
    n = 1;

    while (fabs(term) > TRUNCATION)
    {
      term *= -x / ((double)n);

      sum += term;
      n++;
    }

    printf("%12f %12g %12g %5d\n", x, exp(-x), sum, n - 1);
  }
}

Jak uzyskać lepszą dokładność?

Zauważmy, że $e^{-x} = \frac{1}{e^x}$

Wyznaczenie wartości $e^x$ jest możliwe z większą dokładnością bo kolejne wyrazy ciągu mają ten sam znak.

$$e^{x}=1+x+\frac{x^{2}}{2}+\frac{x^{3}}{6}+\cdots+\frac{x^{n}}{n !}+\cdots = \sum_{n=0}^{\infty} \frac{x^{n}}{n !} $$

exp_recursive_fixed.c

#include <stdio.h>
#include <math.h>

#define TRUNCATION 1e-10
int main()
{
  int loop, n;
  double x, term, sum;

  printf("       x       exp(-x)       series       terms\n");

  for (loop = 0; loop <= 100; loop += 10)
  {
    x = (double)loop;

    sum = 1.0;
    term = 1;
    n = 1;

    while (fabs(term) > TRUNCATION)
    {
      term *= x / ((double)n);

      sum += term;
      n++;
    }

    printf("%12f %12g %12g %5d\n", x, exp(-x), 1.0 / sum, n - 1);
  }
}

Zadanie 2: Funkcja $x - \sin{x}$

Zaimplementuj program wyznaczający wartość funkcji $f(x) = x - \sin{x}$ korzystając z rozwinięcia potęgowego w taki sposób aby zminimalizować utratę precyzji obliczeń oraz wpływ błędów zaokrągleń.

Rozwinięcie funkcji $\sin{x}$ jest postaci: $$\sin x=x-\frac{x^{3}}{3 !}+\frac{x^{5}}{5 !}-\frac{x^{7}}{7 !}+$$

Zaimplementuj obliczenia unikając wyznaczenia wartości silni w sposób bezpośredni.

Wyprowadzenie wzoru rekurencyjnego:

$$x-\sin x=\frac{x^{3}}{3 !}-\frac{x^{5}}{5 !}+\frac{x^{7}}{7 !}-\frac{x^{9}}{9 !}+\ldots = \sum^{N}_{i=1} s_i$$

$$s_{1}=\frac{x^{3}}{3 !}$$

$$s_{2} = -\frac{x^{5}}{5 !} = -\frac{x^{3}}{3 !} \frac{x^{2}}{4 \cdot 5} = -s_{1} \frac{x^{2}}{4 \cdot 5} $$

$$s_{3} =\frac{x^{7}}{7 !} = \frac{x^{5}}{5 !} \frac{x^{2}}{6 \cdot 7}= s_{2} \frac{x^{2}}{6 \cdot 7}$$

$$s_{n}=s_{n-1}(-1)^{n-1} \frac{x^{2}}{2 n(2 n+1)}, \quad n=1, \ldots, N$$ $$s_0 = x$$

Dokładność wyznaczenia funkcji $sin(x)$ zanika dla wartości $x$ odległych od $0$. Postaraj się zniwelować wpływ tego problemu.

Wskazówka: zauważ, że $$ \sin(x) = \sin(x - 2 \pi)$$ oraz $$ \sin(x) = \sin(\pi - x) \quad \text{dla} \quad x \in [\frac{\pi}{2}, \pi]$$

Użyj wartości $\pi$ zdefiniowanej w pliku nagłówkowym biblioteki matematycznej (M_PI) lub przyjmij $\pi \approx 3.14159265358979323846$

Program powinien umożliwić wprowadzenie użytkownikowi wartości rzeczywistej $x$ a następnie wyświetla następujące wartości:

• wartość funkcji $f(x) wyznaczonej z rozwinięcia w szereg Taylora
• wartość dokładną wyznaczoną za pomocą funkcji matematycznych z biblioteki math.h
• wartość błędu bezwzględnego
• wartość błędu względnego (w procentach).

Dodatkowe ćwiczenia

Dodatkowe ćwiczenia nie podlegające ocenie.

Ćw. 1

Napisz program, który wyznacza przybliżoną wartość funkcji $\ln{(1-x)}$ dla $|x| < 1$, korzystając z rozwinięcia potęgowego $$\ln{(1-x)} \approx - \left( x + \frac{x^2}{2} + \frac{x^3}{3} + \frac{x^4}{4} + \ldots \right)$$

Dane wejściowe wprowadzane przez użytkownika:

• wartość $x$
• dokładność obliczeń $\epsilon > 0$

Program wypisuje następujące wartości:

• przybliżona wartość funkcji $\ln{(1-x)}$ z dokładnością $\epsilon$, tj. obliczenia przerywane są gdy kolejny człon rozwinięcia potęgowego spełni $a_i <\epsilon$
• dokładną wartość funkcji wyznaczoną za pomocą funkcji z biblioteki matematycznej
• wartość błędu względnego

Ćw. 2

Napisz program, który doda do siebie $n$ razy liczbę $x$. Program dla podanej przez użytkownika liczby $x > 0$ wypisuje wynik sumowania $\sum_{i=1}^n{x}$ dla $n=10, 100, 1000, 10000, 1000000$. Dla każdej wartości $n$ wypisywana jest też wartość dokładna $n\cdot x$ oraz błąd względny wyniku.

Sprawdź wyniki programu uzyskane dla $x=0.1$ oraz $x=0.25$.