Laboratorium Metod Numerycznych w Inżynierii

Ćwiczenie 1

„Przenoszenie się błędów w obliczeniach numerycznych”

- sprawozdanie z wykonanego ćwiczenia

Autorzy:

Kosiński Tomasz

Leśnik Konrad

Zadanie 1

Listing M-funkcji obliczającej wartość wielomianu bezpośrednio z typowej postaci wielomianu:

function result = postWiel(WSP, vX)

result = 0; % zerowanie wektora wyjściowego

n = length(WSP);

for i = 1: n % obliczenia dla wszystkich współczynników

result = result + WSP(i) * vX^(n-1);

end

print result

Listing M-funkcji obliczającej wartość wielomianu przy pomocy schematu Hornera:

function result = schHorn(WSP, vX)

result = 0;

n = length(WSP) - 1;

for i = 1: n

result = result + WSP(i); % dodanie kolejnego współczynnika (a_n)

result = result * vX; % mnożenie przez wartość argumentu (x)

end

result = result + WSP(n+1); % dodanie ostatniego współczynnika (a₀)

Listing zadanego wielomianu (x-1) ¹⁰:

function result = wiel(X)

result = (X - 1).^1

Listing skryptu wykreślającego wartości otrzymane przy pomocy napisanych

M-funkcji:

WSP = poly(ones(1, 10)); % wektor zawierajacy wspólczynniki wielomianu

vX = linspace(0.95, 1.05); % wektor wartosci x

for i = 1: 100

postacWielomianu(i) = postWiel(WSP, vX(i));

schematHornera(i) = schHorn(WSP, vX(i));

wielomian(i) = wiel(vX(i));

end

plot(vX, postacWielomianu, 'b-', vX, schematHornera, 'r:', vX, wielomian, 'g--') % rysujemy wykres, poniżej opis wykresu

title 'Wartosci otrzymane przy pomocy M-funkcji'

xlabel 'Wartosci argumentow (vX)'; ylabel 'Wartosci wspolczynnikow (WSP)'

legend('Postac wielomianu', 'Schemat Hornera', 'Wielomian')

Wnioski:

Na podstawie otrzymanego wykresu możemy łatwo zaobserwować, że wyniki uzyskane przy pomocy schematu Hornera `lepiej' przybliżają nasz wielomian. Dzieje się tak, gdyż obliczając wartości wielomianu n-tego stopnia bezpośrednio z typowej postaci wykonujemy n*(n+1)/2 mnożeń, natomiast używając schematu Hornera jedynie n mnożeń, co daje znaczne zmniejszenie liczby wykonywanych działań. Choć oba sposoby powinny dać w rezultacie te same wyniki, nie stało się tak w tym przypadku ze względu na ograniczoną reprezentację liczb zmiennoprzecinkowych, a co za tym idzie - metoda o większej liczbie wykonanych działań obarczona jest większym błędem (błąd zaokrągleń).

Zadanie 2

Listing M-funkcji obliczającej przybliżoną wartość pochodnej korzystając ze wzoru na różnicę centralną:

function result = roznicaCentr(x0, h)

result = 0;

result = (cos(x0+h) - cos(x0-h)) / (2*h);

Listing M-funkcji obliczającej przybliżoną wartość pochodnej korzystając ze wzoru na różnicę progresywną:

function result = roznicaProgr(x0, h)

result = 0;

result = (cos(x0+h) - cos(x0)) / h;

Listing M-funkcji obliczającej błąd całkowity dla różnicy centralnej:

function result = bladRc(h)

result = (((h^2)/6)*(sin(h))) / (-sin(h));

Listing M-funkcji obliczającej błąd całkowity dla różnicy progresywnej:

function result = bladRp(h)

result = ((h/2)*(-cos(h))) / (-sin(h));

Listing skryptu wynikowego wykreślającego na jednym wykresie błędy całkowite i względne w skali logarytmicznej:

x0 = pi/15; % wartość początkowa

h = pi/3; % krok

wp1='-sin(x0)'; p1=eval(wp1); % dla ułatwienia obliczamy wartości

wp2='-cos(x0)'; p2=eval(wp2); % pochodnych i przypisujemy je

wp3='sin(x0)'; p3=eval(wp3); % do zmiennych p1, p2, p3

% (dla lepszej czytelności kodu)

for i = 1: 50 % Obliczenia wykonywane dla 50 punktów

H(i) = h*2 ^(-(i-1)); % Wyznaczany krok obliczeń

rc_vals(i) = roznicaCentr(x0,H(i)); % Różnica centralna (przybliżenie pochodnej)

rp_vals(i) = roznicaProgr(x0,H(i)); % Różnica progresywna (j.w.)

rc_blad(i) = ((((H(i))^2)/6)*p3)/p1; % Błąd całkowity z różnicy centralnej

rp_blad(i)=((H(i)/2)*p2)/p1; % Błąd całkowity z różnicy progresywnej

rc_bw(i)=(rc_vals(i)-p1)/p1; % Błąd względny różnicy centralnej

rp_bw(i)=(rp_vals(i)-p1)/p1; % Błąd względny różnicy progresywnej

end

loglog(H, abs(rc_blad), 'rd', H, abs(rp_blad), 'bp', H, abs(rc_bw), 'black--', H, abs(rp_bw), 'm--'); % Wykres w skali logarytmicznej

title 'Wykres bledow calkowitych i wzglednych dla roznicy centralnej i progresywnej'

xlabel 'Krok obliczen'; ylabel 'Wartosci bledow'

legend ('Blad calkowity z roznicy centralnej', 'Blad calkowity z roznicy progresywnej', 'Blad wzgledny roznicy centralnej', 'Blad wzgledny roznicy progresywnej')

Na podstawie otrzymanych wykresów oszacowany przez nas (poprzez rzutowanie punktów przecięcia wykresów błędów całkowitych i względnych dla każdej z metod) optymalny krok obliczeń dla metody progresywnej wynosi ok. 10^-8, natomiast dla różnicy centralnej ok. 10^-5 (wyniki w skali logarytmicznej). Im mniejszy jest krok obliczeń, tym wykresy błędów pokrywają się szybciej; w naszym przypadku lepszym okazał się krok h dla różnicy progresywnej i wynosi on h_OPT = 10^-8

(gdyż błąd względny szybciej `zbiega' do wykresu błędu całkowitego).

---