Nevíte-li si rady s jakýmkoliv matematickým problémem, toto místo je pro vás jako dělané.
Nástěnka
❗22. 8. 2021 (L) Přecházíme zpět na doménu forum.matweb.cz!
❗04.11.2016 (Jel.) Čtete, prosím, před vložení dotazu, děkuji!
❗23.10.2013 (Jel.) Zkuste před zadáním dotazu použít některý z online-nástrojů, konzultovat použití můžete v sekci CAS.
Nejste přihlášen(a). Přihlásit
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.
Offline
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ů.
Offline
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.
Offline
↑ rss:Ahoj, jeden stroj mi vyhodil nasledovné, možno Ti to bude vhod =>
Offline
↑ 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
z rovnice
(viz předpoklady
oné metody), což u kubického polynomu není až tak těžké.
Offline
↑ 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:...
Offline
↑ Eratosthenes: Ahoj, smekám pred tvorcami programov v symbolike .
Offline
↑ 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é.
Offline
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.
Offline
Zdravím, i když nechceš numerickou metodu, tak ti ji nabídnu. Aspoň se podívej jak je jednoduchá:
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 FunctionPř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).
Offline
↑ 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š...
Offline
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.
Offline
↑ 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.
Offline
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é.
Offline
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.
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 FunctionFunkce 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.
Offline
↑ rss:
V Pascalu
Výstup pole stringů (kořenů)
//ř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;Offline