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 błedu 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$

Błąd zaokrąglenia liczb

Liczba_zmiennoprzecinkowa

$$x = \pm m \cdot 2^{c}$$

Błąd względny wynikający ze skonczonej liczby bitów $t$ mantysy $m$

$$\epsilon \le 2^{-t} $$

float $t=23$ bity $ 2^{-t} \approx 1.2 \cdot 10^{-7}$, double $t=52$ bity $ 2^{-t} \approx 2.2 \cdot 10^{-16}$

Błędy operacji arytmetycznych

$$ x_1 \pm x_2 = \frac{x_1 \epsilon \pm x_2 \epsilon}{x_1 \pm x_2} $$

Możliwy znaczny wzrost błędu przy odejmowaniu dużych bliskich sobie liczb

Przykład: miejsca zerowe paraboli

Napisz program znajdujący pierwiastki równania kwadratowego $$y = ax^2 + bx + c$$ dla dowolnych rzeczywistych $a, b, c$. Dla wyznaczonych pierwiastków wypisz wartości funkcji $f(x_1)$ oraz $f(x_2)$ aby zweryfikować poprawność uzyskanego wyniku.

Sprawdź, czy program poprawnie wyznacza miejsca zerowe równania $y = x^2 - 6.433 x + 0.009474$ oraz $y = 0.1 x^2 - 100 x + 0.1$.

Przy wyznaczaniu pierwiastków tradycyjnym podejściem (z wyznaczeniem delty) może dojść do istotnego błędu zaokrągleń gdy delta $\sqrt{\Delta}$ jest bliska $b$. Spróbuj zniwelować wpływ tego błędu wykorzystując Wzory Viète’a

pierwiastki.c

#include <stdio.h>
#include <math.h>
#include <float.h>
 
#define EPS FLT_EPSILON
 
float wielomian(float x, float a, float b, float c)
{
   return x*x*a + x*b + c; 
}
 
int main(void)
{
   float a, b, c, delta, x1, x2 ;
 
   printf("a="); scanf("%f",&a);
   printf("b="); scanf("%f",&b);
   printf("c="); scanf("%f",&c);
 
   printf("Wielomian w(x) = %lg x^2 + %lg x + %lg.\n", a, b, c);
 
   delta = b*b - 4.0*a*c;
 
   if ( delta < 0)
   {
      printf("Brak miejsc zerowych\n");
      return 0;
   }
 
   if( fabs(delta) < EPS )
   {
      x1 = -b/a/2;
      x2 = x1;   
   }
   else  
   {
      delta = sqrt(delta);
      x1 = -(b-delta)/(2*a);   
      x2 = -(b+delta)/(2*a);   
   }
 
   printf("Miejsca zerowe, x1=%lg, x2=%lg \n", x1, x2);
   printf("w(x1)=%lg,\nw(x2)=%lg\n", wielomian(x1, a, b, c), wielomian(x2, a, b , c));
 
   return 0;
}

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

Treść zadania znajduje się w Moodle pod adresem Zadanie 2

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