Matematické Fórum

rss · 19. 06. 2015 09:53

Dobrý deň.

Potrebujem vyjadriť reálne x1 kubickej rovnice pre počítačovú simuláciu s presnosťou aspoň 4 platné čísla, viď obrázok. http://s9.postimg.org/6ztvw9bpp/kubikgravitacia.png . Ak by sa dalo tak aj x2 a x3. Nechcem numerickú metódu.

rss · 19. 06. 2015 09:53

Neviem či je táto otázka dostatočne ťažká a či som ju zaradil správne. :-)

Bati · 19. 06. 2015 10:10

Ahoj,
původní Cardanova metoda najde reálný kořen (pokud tedy jde o rovnici s reálnými koeficienty, jinak nemusí reálný kořen existovat). Hned pod tím v tom odkaze je další hezká metoda pro nalezení všech kořenů.

rss · 23. 06. 2015 20:04

Bati napsal(a):
Ahoj,
původní Cardanova metoda najde reálný kořen (pokud tedy jde o rovnici s reálnými koeficienty, jinak nemusí reálný kořen existovat). Hned pod tím v tom odkaze je další hezká metoda pro nalezení všech kořenů.

Hneď na úvod som použil Cardanovu metódu a vyjadril X všeobecne, Potrebujem vyjadriť len reálne X1.

Bati · 23. 06. 2015 23:54

↑ rss:
Jsem přesvědčen, že tvou otázku jsem zodpověděl (↑ Bati:).

pietro · 24. 06. 2015 11:29

↑ rss:Ahoj, jeden stroj mi vyhodil nasledovné, možno Ti to bude vhod =>

Autor skryl část textu. Zobrazit/skrýt!

Rumburak

↑ rss:

Ahoj.

Soudím, že hledáš nějakou numerickou metodu. Kdysi se mi na obdobnou úlohu (také šlo o hledání reálného kořene
kubické rovnice) osvědčila metoda tečen, která konvergovala velmi rychle. Informace o ní jistě na webu najdeš.
Předpokladem jejího využití je mít dostatečnou představu o průběhu funkce $f$ z rovnice $f(x) = 0$ (viz předpoklady
oné metody), což u kubického polynomu není až tak těžké.

Eratosthenes · 24. 06. 2015 17:12

↑ Rumburak:

jenže ↑ rss: hned na začátku píše "Nechcem numerickú metódu". Škoda. Kdyby ji chtěl, pak takové půlení intervalu je na pět řádků programu a tisícinu vteřiny strojového času. Nechápu proč jsou mu milejší šílenosti typu ↑ pietro:...

pietro · 24. 06. 2015 22:04

↑ Eratosthenes: Ahoj, smekám pred tvorcami programov v symbolike .

rss · 24. 06. 2015 23:03

↑ pietro: Ďakujem, asi na to budem potrebovať viac času, je to to isté čo vyšlo mne?

pietro · 25. 06. 2015 08:22

↑ rss: TenTvoj dáva zaručene vždy reálny koreň ( pre rôzne kombinácie b,d) . Ten "môj reálny" občas vydá ( pri niektorých b,d) aj jeden z komplexných koreňov . Ale keď oba vzorce vydajú reálny koreň, tak sú výsledky u Tvojho a môjho zhodné.

rss · 25. 06. 2015 12:34 — Editoval rss (25. 06. 2015 12:36)

pietro napsal(a):
↑ rss: TenTvoj dáva zaručene vždy reálny koreň ( pre rôzne kombinácie b,d) . Ten "môj reálny" občas vydá ( pri niektorých b,d) aj jeden z komplexných koreňov . Ale keď oba vzorce vydajú reálny koreň, tak sú výsledky u Tvojho a môjho zhodné.

Pri výpočte v OpenOfficeCalc mi vyhadzovalo chybu, záporná odmocnina = komplexné číslo.

mák · 25. 06. 2015 18:13

Zdravím, i když nechceš numerickou metodu, tak ti ji nabídnu. Aspoň se podívej jak je jednoduchá:

Code:

Function CubicFindRealRoot(B#,D#,Z#) As Double
    Dim X As Double 
    X = -9e+99
    While Abs(X-Z)>1e-9
        X = Z
        Z = X-(X^3+B*X^2+D)/(3*X^2+2*B*X)
    Wend
    CubicFindRealRoot = Z
End Function

Předáš koeficienty B a D. Proměnná Z je odhad, čím přesněji jej zadáš, tím se méně cyklů se použije. Vrací nejbližší reálný kořen s přesností na 9 desetinných míst (porovnání While).

Eratosthenes · 25. 06. 2015 20:02

↑ mák:

Dva dotazy:

1) Pokud jsem stačil zkontrolovat, funkce řeší rovnici x^3+bx^2+d=0. Co lineární člen? Obecně je přeci x^3+bx^2+cx+d=0...

2) Jak víš, že to zkonverguje? Zkus dosadit b=2, d=5, z=0.1. Myslím, že se kořene nedobereš...

mák · 25. 06. 2015 22:32

1) Funkce řeší rovnici x^3+B*x^2+D podle zadání uvedeném v prvním příspěvku. Pro obecnou rovnici by byla funkce o něco složitější. Chtěl jsem, aby se výpočet zbytečně nekomplikoval.

2) Pro uvedený příklad b=2, d=5, z=0.1 funkce konverguje k číslu -2,69064744802861, Výsledek (reálný kořen) je podle wxMaximy x=-2.69064744802863. Takže zrovna tady to vyšlo. Je možné, že pro některá čísla to nebude konvergovat. Mám však zkušenost, že Newtonova metoda konverguje většinou velmi dobře a rychle (postačuje pouze několik cyklů) zejména, pokud počáteční nástřel je blízko.

Eratosthenes

↑ mák:

koukám, že k matematice přistupuješ jako k loterii - když při Člověče nezlob se nechceš jedničku, taky hod většinou vyjde. Newtonova metoda konverguje sice rychle, ale podmínky konvergence jsou dost složité. Ten první příklad jsem napsal špatně - zkus rovnici

x^3+5*x^2+3=0

s "počátečním nástřelem" x=-0.5. A to už vůbec nemluvím o "nástřelu" x=0. Takže do programu bych něco takového určitě nedával.

mák · 26. 06. 2015 00:31

Ten tvůj příklad:

Calc (macro): -5,11467914980952
wxMaxima (reálný kořen): x=-5.114679149809521

Takže to zase vyšlo...

K té loterii:
Pokud musím numericky řešit nějaký výpočet, tak vím v jakém rozmezí hodnot budou vstupní a výstupní hodnoty. Takže vytvořím jednoduchý náhradní výpočet, který mi zaručí, že nástřel se bude lišit od skutečnosti pouze minimálně a Newtonovu metodu použiji pouze pro zpřesnění na požadovaný počet míst. Takto to řeším já a zatím mi to vychází.

Proto jsem řekl, že záleží na počátečním odhadu. To bude muset ale vyřešit někdo jiný.

A nebo zkusit jiný způsob výpočtu.
Já jsem jenom chtěl ukázat, že makro je jednoduché.

mák · 26. 06. 2015 04:26 — Editoval mák (26. 06. 2015 04:27)

Ještě se k tomu vracím, nedalo mi to a napsal jsem přímou metodu výpočtu všech třech kořenů včetně komplexních.

Code:

Function acos(X#) As Double 
    acos = Atn(-X / Sqr(-X * X + 1)) + 2 * Atn(1) 
End Function

Function QBRT(X As Double) As Double
    QBRT = Abs(X) ^ (1 / 3) * Sgn(X)
End Function

Function CubicRoots(A#, B#, C#, D#) As Variant
    Dim F#,G#,H#,I#,J#,K#,L#,M#,N#,O#,P#,Q#,R#,S#,T#,U#,W#
    Dim Z(2,1) As Variant, Y As Integer 
    For Y=0 To 2
        Z(Y,0) = ""
        Z(Y,1) = ""
    Next Y     
    F = C/A-(B/A)^2/3
    G = (27*A^2*D-9*A*B*C+2*B^3)/(27*A^3)
    H = G^2/4+F^3/27
    If H<0 Then ' Tři různé reálné kořeny
        I = sqr(G^2/4-H)
        J = QBRT(I)
        K = acos(-G/(2*I))
        M = cos(K/3)
        N = sqr(3)*sin(K/3)
        P = -B/(3*A)
        Z(0,0) = 2*J*M+P
        Z(1,0) = P-J*(M+N)
        Z(2,0) = P-J*(M-N)
    ElseIf H>0 Then ' Jeden reálný kořen a dva komplexní
        R = -G/2
        T = sqr(H)
        S = QBRT(R+T)
        U = QBRT(R-T)
        L = S+U
        P = -B/(3*A)
        Z(0,0) = +L  +P
        Z(1,0) = -L/2+P
        Z(2,0) = Z(1,0)
        Z(1,1) = +(S-U)*sqr(3)/2
        Z(2,1) = -Z(1,1)
    Else ' Tři identické reálné kořeny
        W = QBRT(-(D/A))
        Z(0,0) = W
        Z(1,0) = W
        Z(2,0) = W    
    EndIf
    CubicRoots = Z        
End Function

Funkce vrací pole 3x2 (levý sloupeček reálná hodnota, pravý imaginární) tří kořenů. Funkce se zadává maticově, takže klasicky trojhmatem [Control]+[Shift]+[Return].
Počítá obecnou rovnici A*X^3+B*X^2+C*X+D = 0.

mák · 26. 06. 2015 11:20

Demonstrační program ke stažení zde.

Honzc · 02. 07. 2015 13:31

↑ rss:
V Pascalu
Výstup pole stringů (kořenů)

Code:

//řeší rovnici: ax^3+bx^2+cx+d=0
//do uses přidat knihovnu Math

type 
  TAlgr=array of string;

//pres - zadá se na kolik platných číslic chceme výsledek (číslo, např.8)

function Algr3(var a,b,c,d: Extended; pres: Integer): TAlgr;
var p,q,dis,y1,y2,y3,yr,yi,x1,x2,x3,fi,r1: Extended;
  rc,ic1,ic2,sy,sx,pr: string;
begin
  SetLength(Result,3);
  pr := '%.'+IntToStr(pres)+'g';
  if a=0 then
  begin
    Result[0] := 'Rovnice není 3.stupně';
    Result[1] := '';
    Result[2] := '';
    MessageBeep($FFFFFFFF);
    Exit;
  end;
    p := (3*a*c-b*b)/(3*a*a);
  if abs(p)<1E-14 then
    p := 0;
  q := (2*b*b*b)/(27*a*a*a)-(b*c)/(3*a*a)+(d/a);
  if abs(q)<1E-14 then
    q := 0;
  dis := -Power(p/3,3)-Power(q/2,2);
  if abs(dis)<1E-15 then
    dis := 0;
  if q>=0 then
    r1 := Sqrt(abs(p/3))
  else
    r1 := -Sqrt(abs(p/3));
  if p<0 then
  begin
    if dis>=0 then
    begin
      try
        if abs(q/(2*Power(r1,3))-1)<1E-14 then
          fi := 0
        else
          fi := ArcCos(q/(2*Power(r1,3)));
        y1 := -2*r1*Cos(fi/3);
        if abs(y1)<1E-14 then
          y1 := 0;
        y2 := 2*r1*Cos(Pi/3-fi/3);
        if abs(y2)<1E-14 then
          y2 := 0;
        y3 := 2*r1*Cos(Pi/3+fi/3);
        if abs(y3)<1E-14 then
          y3 := 0;
        sx := Format('%.15g', [b/(3*a)]);
        sy := Format('%.15g', [y1]);
        x1 := StrToFloat(sy)-StrToFloat(sx);
        if abs(x1)<1E-14 then
          x1 := 0;
        sy := Format('%.15g', [y2]);
        x2 := StrToFloat(sy)-StrToFloat(sx);
        if abs(x2)<1E-14 then
          x2 := 0;
        sy := Format('%.15g', [y3]);
        x3 := StrToFloat(sy)-StrToFloat(sx);
        if abs(x3)<1E-14 then
          x3 := 0;
        Result[0] := Format(pr, [x1]);
        Result[1] := Format(pr, [x2]);
        Result[2] := Format(pr, [x3]);
      except
        on EInvalidOp do
        begin
          Result[0] := 'Přetečená paměť!';
          Result[1] := '';
          Result[2] := '';
        end;
      end;
    end
    else
    begin
      try
        if abs(q/(2*Power(r1,3))-1)<1E-14 then
          fi := 0
        else
          fi := ArcCosh(q/(2*Power(r1,3)));
        y1 := -2*r1*Cosh(fi/3);
        if abs(y1)<1E-14 then
          y1 := 0;
        yr := r1*Cosh(fi/3);
        if abs(yr)<1E-14 then
          yr := 0;
        yi := Sqrt(3)*r1*Sinh(fi/3);
        if abs(yi)<1E-14 then
          yi := 0;
        sx := Format('%.15g', [b/(3*a)]);
        sy := Format('%.15g', [y1]);
        x1 := StrToFloat(sy)-StrToFloat(sx);
        if abs(x1)<1E-14 then
          x1 := 0;
        sy := Format('%.15g', [yr]);
        x2 := StrToFloat(sy)-StrToFloat(sx);
        if abs(x2)<1E-14 then
          x2 := 0;
        if (x2=0) then
        begin
          if (yi=0) then
          begin
            rc := Format(pr, [x2]);
            ic1 := '';
            ic2 := '';
          end
          else
          begin
            rc := '';
            ic1 := '- i '+Format(pr, [abs(yi)]);
            ic2 := 'i '+Format(pr, [abs(yi)]);
          end;
        end
        else
        begin
          rc := Format(pr, [x2]);
          if (yi=0) then
          begin
            ic1 := '';
            ic2 := '';
          end
          else
          begin
            ic1 := ' - i '+Format(pr, [abs(yi)]);
            ic2 := ' + i '+Format(pr, [abs(yi)]);
          end;
        end;
        Result[0] := Format(pr, [x1]);
        Result[1] := rc+ic1;
        Result[2] := rc+ic2;
      except
        on EInvalidOp do
        begin
          Result[0] := 'Přetečená paměť!';
          Result[1] := '';
          Result[2] := '';
        end;
      end;
    end;
  end
  else
  begin
    if p<>0 then
    begin
      try
        if abs(q/(2*Power(r1,3)))<1E-15 then
          fi := 0
        else
          fi := ArcSinh(q/(2*Power(r1,3)));
        y1 := -2*r1*Sinh(fi/3);
        if abs(y1)<1E-14 then
          y1 := 0;
        yr := r1*Sinh(fi/3);
        if abs(yr)<1E-14 then
          yr := 0;
        yi := Sqrt(3)*r1*Cosh(fi/3);
        if abs(yi)<1E-14 then
          yi := 0;
        sx := Format('%.15g', [b/(3*a)]);
        sy := Format('%.15g', [y1]);
        x1 := StrToFloat(sy)-StrToFloat(sx);
        if abs(x1)<1E-14 then
          x1 := 0;
        sy := Format('%.15g', [yr]);
        x2 := StrToFloat(sy)-StrToFloat(sx);
        if abs(x2)<1E-14 then
          x2 := 0;
        if (x2=0) then
        begin
          if (yi=0) then
          begin
            rc := Format(pr, [x2]);
            ic1 := '';
            ic2 := '';
          end
          else
          begin
            rc := '';
            ic1 := '- i '+Format(pr, [abs(yi)]);
            ic2 := 'i '+Format(pr, [abs(yi)]);
          end;
        end
        else
        begin
          rc := Format(pr, [x2]);
          if (yi=0) then
          begin
            ic1 := '';
            ic2 := '';
          end
          else
          begin
            ic1 := ' - i '+Format(pr, [abs(yi)]);
            ic2 := ' + i '+Format(pr, [abs(yi)]);
          end;
        end;
        Result[0] := Format(pr, [x1]);
        Result[1] := rc+ic1;
        Result[2] := rc+ic2;
      except
        on EInvalidOp do
        begin
          Result[0] := 'Přetečená paměť!';
          Result[1] := '';
          Result[2] := '';
        end;
      end;
    end
    else
    begin
      if ((a>0)and(d>0))or((a<0)and(d<0)) then
        r1 := -1
      else
        r1 := 1;
      y1 := r1*Power(abs(q),1/3);
      if abs(y1)<1E-14 then
        y1 := 0;
      x1 := y1-b/(3*a);
      if (q<>0) then
      begin
        yr := -y1/2-b/(3*a);
        yi := y1*Sqrt(3)/2;
        ic1 := ' - i '+Format(pr, [abs(yi)]);
        ic2 := ' + i '+Format(pr, [abs(yi)]);
        rc := Format(pr, [yr]);
        Result[0] := Format(pr, [x1]);
        Result[1] := rc+ic1;
        Result[2] := rc+ic2;
      end
      else
      begin
        Result[0] := Format(pr, [x1]);
        Result[1] := Format(pr, [x1]);
        Result[2] := Format(pr, [x1]);
      end;
    end;
  end;
end;

Matematické Fórum

#1 19. 06. 2015 09:53

kubická rovnica, x1 pre výpočet počítačom

#2 19. 06. 2015 09:53

Re: kubická rovnica, x1 pre výpočet počítačom

#3 19. 06. 2015 10:10

Re: kubická rovnica, x1 pre výpočet počítačom

#4 23. 06. 2015 20:04

Re: kubická rovnica, x1 pre výpočet počítačom

Bati napsal(a):

#5 23. 06. 2015 23:54

Re: kubická rovnica, x1 pre výpočet počítačom

#6 24. 06. 2015 11:29

Re: kubická rovnica, x1 pre výpočet počítačom

#7 24. 06. 2015 16:46 — Editoval Rumburak (24. 06. 2015 16:58)

Re: kubická rovnica, x1 pre výpočet počítačom

#8 24. 06. 2015 17:12

Re: kubická rovnica, x1 pre výpočet počítačom

#9 24. 06. 2015 22:04

Re: kubická rovnica, x1 pre výpočet počítačom

#10 24. 06. 2015 23:03

Re: kubická rovnica, x1 pre výpočet počítačom

#11 25. 06. 2015 08:22

Re: kubická rovnica, x1 pre výpočet počítačom

#12 25. 06. 2015 12:34 — Editoval rss (25. 06. 2015 12:36)

Re: kubická rovnica, x1 pre výpočet počítačom

pietro napsal(a):

#13 25. 06. 2015 18:13

Re: kubická rovnica, x1 pre výpočet počítačom

Code:

#14 25. 06. 2015 20:02

Re: kubická rovnica, x1 pre výpočet počítačom

#15 25. 06. 2015 22:32

Re: kubická rovnica, x1 pre výpočet počítačom

#16 25. 06. 2015 23:21 — Editoval Eratosthenes (26. 06. 2015 00:00)

Re: kubická rovnica, x1 pre výpočet počítačom

#17 26. 06. 2015 00:31

Re: kubická rovnica, x1 pre výpočet počítačom

#18 26. 06. 2015 04:26 — Editoval mák (26. 06. 2015 04:27)

Re: kubická rovnica, x1 pre výpočet počítačom

Code:

#19 26. 06. 2015 11:20

Re: kubická rovnica, x1 pre výpočet počítačom

#20 02. 07. 2015 13:31

Re: kubická rovnica, x1 pre výpočet počítačom

Code:

Zápatí