Wyznacza szybko i dokładnie rozwiązania rzeczywiste, z ich krotnością

Mateusz Kowalski

Wyprowadzenie wzorów matematycznych na równanie 3 stopnia.

Krok po kroku W poniższym filmie po polsku.

Spis Treści nagrania
3:50 Sprowadzenie zagadnienie do równania 3 stopnia bez składnika bez niewiadomej w kwadracie.
7:40 wprowadzenie pomocniczych niewiadomych u i v
10:35 Sprowadzenie do równania trikwadratowego
14:55 pierwiastki z jedynki (dygresja)
21:35 własności pierwiastków z 1
22:00 pierwiastki zespolone z liczby rzeczywistej
24:40 powrót do zadania
24:40 wyznaczenie u i v
33:52 końcowe wzory dla delty większej lub równej 0
35:24 obserwacje wyników
37:37 Gdy delta równa się 0 pierwiastki podwójne i potrójne
39:40 Gdy delta mniejsza od 0
41:00 pierwiastkowanie liczby zespolonej
41:52 płaszczyzna zespolona
45:55 pierwiastki zespolone wzory
47:50 wyznaczenie u i v
56:00 końcowe wzory dla delty mniejszej od 0
57:15 sin 10 stopni dokładna wartość
57:30 wyprowadzenie wzoru na sinus potrójnego kąta
59:00 Przykład użycia wzorów

Omówienie Implementacji algorytmu oraz rozwiązań potencjalnych problemów i poprawą precyzji obliczeń

Krok po kroku W poniższym filmie po polsku.

Spis Treści nagrania
2:15 Plan nagrania
4:42 Wzory rozwiązujące równanie 3 stopnia
7:35 Ciekawy problem
12:07 Kosinus z 1/3 kąta cos(arcocso(x)/3)
15:40 Cenna uwaga
25:14 Czy naprawdę jest to problem
31:32 Algorytm rozwiązujący
34:42 Pisanie algorytmu
44:17 Potencjalne problemy
45:40 Wyróżnik równania 3 stopnia
46:22 Dalsze problemy
52:35 Wpływ q
54:52 Dobry wzór na fi
57:25 Algorytm drugie podejście
1:04:22 Dalsze udoskonalenia
1:06:30 Dlaczego nie dla metody p przez q
1:07:45 Równanie zwrotne
1:08:30 Poprawa błędów obliczeniowych
1:13:50 Poprawa funkcji sqrt3
1:17:30 Trzecia wersja algorytmu
1:51:52 Testowanie algorytmu
1:53:02 Poprawa dokładności rozwiązań równania

Kody Algorytmu napisane, podczas filmu w środowisku obliczeniowym Matlab / Octave

Funkcja obliczająca pierwiastek arytmetyczny 3 stopnia, również dla liczb ujemnych, z poprawą błędów numerycznych.

function x = sqrt3(a)
  if a >= 0 
    x = a^(1/3);
  else
    x = -(-a)^(1/3);
  end
  if a~=0 
    x = x - (x^3-a)/(3*x^2);
  end
end 

Funkcja rozwiązująca równanie 3 stopnia ax3 + bx2 + cx + d = 0 o współczynnikach rzeczywistych, również dla a=0

function [x, krot] = row3stv3(a,b,c,d)
  if a == 0
    if b==0
      if c ==0
        if d == 0
          error('równanie nic nie wnosi');
        else
          error('równanie sprzeczne');
        end
      else
        x = [-d/c NaN NaN];
        krot = [1 0 0];
        return;
      end
    else
      del = c^2 - 4*b*d;
      if del > 0
        x = [(-c-sqrt(del))/(2*b) (-c+sqrt(del))/(2*b) NaN];
        krot = [1 1 0];
        return;        
      elseif del < 0
        x = [NaN NaN NaN] ;
        krot = [0 0 0 ];
        return;
      else
        x = [-c/(2*b) NaN NaN];
        krot = [2 0 0 ];
        return;
      end
    end
  end
  
  % od tego miejsca wiadomo że a~=0
  if d ==0
    del = b^2 - 4*a*c
    if del > 0
      if -b-sqrt(del) == 0
        x = [0 (-b+sqrt(del))/(2*a) NaN ];
        krot = [2 1 0];
        return;
      elseif -b+sqrt(del) == 0
        x = [0 (-b-sqrt(del))/(2*a) NaN];
        krot = [2 1 0];
        return;
      else
        x = [0 (-b-sqrt(del))/(2*a) (-b+sqrt(del))/(2*a)];
        krot = [1 1 1];
        return;
      end
    elseif del <0
      x = [0 NaN NaN]
      krot = [1 0 0 ];
      return;
    else
      if b == 0
        x  = [0 NaN NaN ];
        krot = [3 0 0];
        return;
      else
         x  = [0 (-b)/(2*a) NaN ];
        krot = [1 2 0];
        return;
      end
    end
  end
  
  % wiadomo już że a~= 0 oraz d~=0
  if a == d && b == c
    del = b^2 - 2*a*b -3*a^2;
    if del > 0 
      if ((-b+a) + sqrt(del)) == -(2*a)
        x = [-1 ((-b+a) - sqrt(del))/(2*a) NaN];
        krot = [2 1 0];
        return;
      elseif ((-b+a) - sqrt(del)) == -(2*a)
        x = [-1 ((-b+a) + sqrt(del))/(2*a) NaN];
        krot = [2 1 0];
        return;
      else
        x = [-1 ((-b+a) + sqrt(del))/(2*a) ((-b+a) - sqrt(del))/(2*a) ];
        krot = [1 1 1];
        return;
     end   
    elseif  del < 0
      x = [-1 NaN NaN];
      krot = [1 0 0];
      return;
    else 
      if (b-a) == 2*a
        x = [-1 NaN NaN];
        krot = [3 0 0];
        return;
      else
        x = [-1 (-b+a)/(2*a) NaN];
        krot = [1 2 0];
        return;
      end
    end
  end
  
  % od tego miejsca już a~= 0 oraz d~=0 oraz (a~=d || b~=c)
  jak_p = -b^2 + 3*a*c;
  jak_q =  2*b^3 - 9*a*b*c + 27*a^2*d;
  p =  (c/a) - (b^2)/(3*a^2);
  q =  (2*b^3)/(27*a^3) + d/a -(b*c)/(3*a^2);
  
  z_y_na_x = -b/(3*a);
  if jak_p == 0
    if jak_q ==0
      x = [z_y_na_x NaN NaN];
      krot = [3 0 0];
      return;
    else 
      x = [sqrt3(-q)+z_y_na_x NaN NaN];
      krot = [1 0 0 ];
	  return;
    end
  else
    if jak_q ==0
      if p > 0
        x =[z_y_na_x NaN NaN];
        krot =[1 0 0 ];
        return;
      else 
        x =[0 sqrt(-p) -sqrt(-p)]+z_y_na_x;
        krot = [1 1 1 ];
        return;
      end
  end
  
  % od tego miejsca już a~= 0 oraz d~=0 oraz (a~=d || b~=c) oraz p~=0 oraz q~=0
  % przypadek ogólny 
  
  eps = 0.01;
  del = q^2 + 4*p^3/27;
  znak_del = 27*a^2*d^2 - 18*a*b*c*d + 4*a*c^3 + 4*b^3*d - b^2*c^2;
  if znak_del > 0
    y(1,1) = sqrt3( (-q-sqrt(del))/2 ) + sqrt3( (-q+sqrt(del))/2 );
    krot = [1 0 0];
  elseif znak_del < 0
    if q > eps
      fi = atan(sqrt(-del)/q) + pi;
    elseif q >= -eps && q <= eps
      fi = acos( ((3*q)/(2*p))*sqrt(3/-p) );
    else
      fi = atan(sqrt(-del)/q) ;
    end
    k = [0 1 2];
    y = 2*sqrt(-3*p)*cos( (fi + 2*k*pi)/3 ) /3;
    krot = [1 1 1];
  else
    y(1,1) = 2*sqrt3(-q/2);
    y(1,2) = sqrt3(q/2);
    krot = [1 2 0];
  end  
  x = y -b/(3*a);
  % poprawka błedów numerycznych
  x = x -  (a*x.^3 + b*x.^2 + c*x + d)./(3*a*x.^2 + 2*b*x + c);
end

Bardzo jestem ciekaw Twojej opinii o tych materiałach. Zostaw komentarz poniżej.

Pozdrawiam
Mateusz Kowalski

0 komentarzy