częste błędy w pracy domowej po zajęciach 5

2020.05.26 — dotyczy: 2019/2020 lato, rw

% Poniżej, aż do odwołania jest przykład błędnego rozwiązania, 
% które z pozoru wygląda sensownie. W dalszej części komentarzy 
% zademonstruję w czym jest błąd.  
% Początek przesłanej pracy
  clear
  p = [2 1 -3 7 -9];
  a = -3;
  b = 0;
  %tolerancja
  t = 1e-4;
  while abs(a-b) > t
    c = (a+b)/2;
    if abs(polyval(p,c)) <= t;
      break
    elseif polyval(p,c)*polyval(p,a) < 0;
      b = c;
    else 
      a = c;
    end
  end
  disp((a+b)/2) 
% koniec oryginalnie przesłanej pracy do sprawdzenia
% drobnym może nie błędem, ale rzeczą niepotrzebną jest liczenie 
% i wyświetlanie (a+b)/2 (linia 18 org. rozwiązania) ponieważ 
% wcale nie musi być to bliższa wartość niż wcześniej policzone c. 
% Także disp(c) jest równie dokładne. Ponadto disp(a+b)/2 
% jest/może być liczone na podstawie nowych wartości a lub 
% b (linia 13 i 15 org. rozwiązania). To jest kwestia do 
% przemyślenia dla wszystkich. Nie traktuje tego jako błąd, ale 
% jest to zdecydowanie zbędne.
% Ale nie o tym. Przejdźmy do ciekawego problemu, którego nie 
% przewidział autor programu (co ciekawe, a raczej zastanawiające, 
% że dokładnie ten sam problem, wystąpił u 30(!) osób. Przypadek?
% Jeżeli ktoś uruchomi ten przykład (albo cały ten skrypt) to 
% zauważy, że program znajduje pierwiastek z zadaną dokładnością, 
% ale działa to tylko dla tego konkretnego przykładu i nie 
% realizuje dobrze bisekcji. W zadaniu jasno było powiedziane, że 
% dokładność pierwiastka powinna być określona i przyjmijmy tu 
% rzeczywiście 1e-4.
% W oryginalnym programie test w linii 10 i wyjście z pętli 
% w linii 11 jeżeli ten test jest prawdziwy jest błędny, ponieważ 
% on sprawdza czy wartość funkcji jest jest bliska zeru dla 
% zadanego argumentu c. Ale jeżeli funkcja jest płaska w okolicach 
% szukanego pierwiastka to ten algorytm po prostu kompletnie 
% zawodzi.
% Kontr-przykład, przejaskrawiony na potrzebę prezentacji problemu.
% Weźmy wielomian W(x) = 0.00001x^2 - 0.0001
% Każdy łatwo w pamięci obliczy że sprowadza się to
% do równania x^2 = 1, a więc pierwiastki to -1 i 1.
% zastosujmy algorytm w granicach od -3 do 0. Powinniśmy
% znaleźć c = -1
% sprawdźmy:
clear
p = [1e-5 0 -1e-5];
a = -3;
b = 0;
%tolerancja
t = 1e-4;
while abs(a-b) > t
  c = (a+b)/2;
  if abs(polyval(p,c)) <= t;
    break
  elseif polyval(p,c)*polyval(p,a) < 0;
    b = c;
  else 
    a = c;
  end
end
c
% c wyszło równe -1.5, czyli błąd tego pierwiastka to
% 0.5, znacznie powyżej
% A teraz zmodyfikujmy nasz program do poprawnego rozwiązania
clear
p = [1e-5 0 -1e-5];
a = -3;
b = 0;
%tolerancja
t = 1e-4;
a = -3;
b = 0;
while abs(a-b) > t
  c = (a+b)/2;
  if polyval(p,c)*polyval(p,a) < 0;
    b = c;
  else 
    a = c;
  end
end
c
% i teraz c jest równe -1 i takiej wartości szukaliśmy
[powrót]