PolyRoot
© 2011 by Andrzej Okuniewski
// zaokrąglanie do zadanej ilości cyfr znaczących function RoundSigDigs($number, $sigdigs) { if($number < 0) { $number = $number * -1; $multiplier = -1; } else { $multiplier = 1; } while ($number < 0.1) { $number *= 10; $multiplier /= 10; } while ($number >= 1) { $number /= 10; $multiplier *= 10; } return round($number, $sigdigs) * $multiplier; } // znak liczby function sign($number) { if($number > 0) { return 1; } else if($number < 0) { return -1; } else if($number == 0) { return 0; } else { return null; } } // ujemna reszta z dzielenia wielomianu p przez q function irem($po, $qo) { $np = count($po) - 1; $nq = count($qo) - 1; $nd = $np - $nq; for($i = 0; $i <= $np; $i++) { $d[$i] = 0; } while($np >= $nq) { $d[$nd] = $po[$np] / $qo[$nq]; for($i = 0; $i <= $nq; $i++) { $po[$np - $i] -= $d[$nd] * $qo[$nq - $i]; } $np--; } $return = null; for($i = $np; $i >= 0 ; $i--) { if(($po[$i] != 0) || (count($return) > 0)) { $return[$i] = -$po[$i]; } } return $return; } //wartość wielomianu w punkcie x function val($po, $x) { $val = 0; for($i = 0; $i < count($po); $i++) { $val += $po[$i]*pow($x, $i); } return $val; } //liczba zmian znaków w ciągu Sturma function nu($po, $x) { echo "x = $x: { "; $nu = 0; for($i = 0; $i < (count($po)-1); $i++) { $sign = sign(val($po[$i], $x)); if($i > 0) { if($last_sign != 0 && $sign != 0) { if($last_sign != $sign) { $nu++; } } } if($sign != 0) { $last_sign = $sign; } echo $sign." "; } echo " }, ν($x) = $nu
"; return $nu; } $var_set = 0; $a = null; // sprawdzanie podanych zmiennych i tworzenie tablicy współczynników if(isset($_GET["a"])) { if($_GET["a"] != "") { $var_set = 1; $a = trim($_GET["a"]); while(strpos($a, " ") > 0) { $a = str_replace(" ", " ", $a); } $a = array_reverse(explode(" ", trim($a))); while($a[count($a)-1] == 0) { unset($a[count($a)-1]); } } } if($var_set > 0) { // czy zmienne są podane $n = count($a) - 1; // stopień wielomianu echo "Pokaż pełne rozwiązanie
Poszukiwane są pierwiastki poniższego wielomianu $n stopnia:
"; echo "p(x) = "; for($i = $n; $i > 0; $i--) { echo $a[$i]." x"; if($i > 1) { echo "$i"; } echo " + "; } echo $a[0]."
"; // pierwsza pochodna for ($i = $n-1; $i >= 0; $i--) { $ap[$i] = $a[$i+1] * ($i+1); } // druga pochodna for ($i = $n-2; $i >= 0; $i--) { $app[$i] = $ap[$i+1] * ($i+1); } echo "Obliczamy kolejne wyrazy ciągu Sturma pi dla wielomianu p:
"; echo "p0(x) = p(x) =
= ";
$po[0] = $a;
for($i = count($po[0]) - 1; $i > 0; $i--) {
echo $po[0][$i]." x";
if($i > 1) {
echo "$i";
}
echo " + ";
}
echo $po[0][0]."
p1(x) = p′(x) =
= ";
$po[1] = $ap;
for($i = count($po[1]) - 1; $i > 0; $i--) {
echo $po[1][$i]." x";
if($i > 1) {
echo "$i";
}
echo " + ";
}
echo $po[1][0]."
p$iter(x) = −rem(p".($iter-2).", p".($iter-1).") =
= ";
$po[$iter] = irem($po[$iter-2], $po[$iter-1]);
for($i = count($po[$iter]) - 1; $i > 0; $i--) {
echo $po[$iter][$i]." x";
if($i > 1) {
echo "$i";
}
echo " + ";
}
echo $po[$iter][0];
if(count($po[$iter]) == 0) {
$is_rem = false;
$pc = $iter;
echo "0";
} else if($po[$iter][0] == 0) {
$is_rem = false;
$pc = $iter;
}
echo "
gdzie p′ to pierwsza pochodna " ."p, a rem(p, q) to reszta z dzielenia " ."wielomianu p przez wielomian q.
"; echo "Obliczamy znaki wartości kolejnych wyrazów ciągu Struma w −∞ i +∞ oraz ilość zmian tych znaków (ν):
"; echo "x = −∞: { "; $num = 0; for($i = 0; $i < $pc; $i++) { $sign = sign($po[$i][count($po[$i])-1]); if((count($po[$i])-1) & 1) { $sign = -$sign; } if($i > 0) { if($last_sign != 0 && $sign != 0) { if($last_sign != $sign) { $num++; } } } $last_sign = $sign; echo $last_sign." "; } echo " }, ν(−∞) = $num
"; echo "x = +∞: { "; $nup = 0; for($i = 0; $i < $pc; $i++) { $sign = sign($po[$i][count($po[$i])-1]); if($i > 0) { if($last_sign != 0 && $sign != 0) { if($last_sign != $sign) { $nup++; } } } $last_sign = $sign; echo $last_sign." "; } echo " }, ν(+∞) = $nup
"; $zc = $num-$nup; // ilość pierwiastków if($zc == 0) { echo "Ponieważ ν(−∞) − " ."ν(+∞) = 0, to wielomian p nie ma " ."pierwiastków rzeczywistych, co kończy algorytm.
" ."Końcowy wynik
Wielomian nie ma pierwiastków rzeczywistych.
"; } else { echo "Na podstawie twierdzenia Sturma wiadomo, że całkowita liczba " ."rzeczywistych pierwiastków wielomianu p wynosi " ."ν(−∞) − ν(+∞) = $zc. "; echo "
Przedział, w którym znajdują się wszystkie rzeczywiste pierwiastki " ."wielomianu dany jest następującą nierównością Cauchy'ego:
"; echo "|x| ≤ 1 + max " ."0 ≤ i ≤ n−1 " ."|ai| / |an|
"; $M = 0; for($i = 0; $i < $n; $i++) { $temp = abs($a[$i]/$a[$n]); if($temp > $M) { $M = $temp; } } $M++; echo "Więc:
"; echo "x ∈ [−$M, $M]
"; echo "Obliczenie wartości ν na końcach tego przedziału:
"; $num = nu($po, -$M); $nup = nu($po, $M); if($num-$nup == $zc) { echo "pozwala jednoznacznie stwierdzić, że w przedziale tym mieszczą się wszystkie ($zc) pierwiastki rzeczywiste wielomianu p.
"; } else { echo "pozwala stwierdzić, że niektóre pierwiastki znajdują się poza tym przedziłem. Jest to błąd w obliczeniach!
"; } if($zc == 1) { echo "Ponieważ w przedziale tym jest tylko jeden pierwiastek " ."można od razu przejść do jego udokładniania.
"; $boundA[1] = -$M; $boundB[1] = $M; } else { echo "Ponieważ w przedziale tym jest więcej niż jeden pierwiastek, to " ."trzeba wyznaczyć $zc mniejsze podprzedziały, w których będzie " ."tylko jeden pierwiastek (tzw. przedziały odosobnienia). " ."W tym celu będziemy połowić przedział [−$M, $M] aż do " ."uzyskania wszystkich przedziałów odosobnienia.
"; $cut = 0; while($odos < $zc) { $odos = 0; $cut++; $W = 2*$M*pow(2, -$cut); $nupr = $num; $boundApr = -$M; $boundBpr = $M; $nupr = nu($po, -$M); for($i = 1; $i <= pow(2, $cut); $i++) { $nu = nu($po, -$M + $i*$W); $boundBpr = -$M + $i*$W; if(($nupr-$nu) == 1) { $odos++; $boundA[$odos] = $boundApr; $boundB[$odos] = $boundBpr; echo "
*
"; } $boundApr = $boundBpr; $nupr = $nu; } echo "Przy podziale przedziału [−$M, $M] na ".(pow(2, $cut))." części udało się wyznaczyć $odos z $zc przedziałów odosobnienia"; if($odos > 0) { echo " – przedziały te oznaczono gwiazdką.
"; } else { echo "."; } echo ""; } echo "
Na podstawie powyższych rozważań udało się stwierdzić, " ."że pierwiastki wielomianu p leżą w następujących " ."przedziałach:
"; for($i = 1; $i <= $zc; $i++) { echo "x$i" ." ∈ [$boundA[$i], $boundB[$i]]
"; } } // koniec wyszukiwania oddziałów odosobnienia // rozpoczynamy udokładnianie pierwiastków $xc = 0; // ilość udokładnionych pierwiastków for($j = 1; $j <= $zc; $j++) { //j - numer pierwiastka echo "Pierwiastek x$j
"; echo "Pierwiastek x$j leży w przedziale [$boundA[$j], $boundB[$j]]. Do metody Newtona obieramy środek tego przedziału:
"; echo "ξ0 = ".(($boundA[$j]+$boundB[$j])/2)."
"; $xi = (($boundA[$j]+$boundB[$j])/2); $halt = 100; $epsilon = 100; echo "Dopóki εn > 10−8 metodą Newtona obliczamy:
"; echo "n | ξn = ξn−1 − p(ξn−1) / p′(ξn−1 ) | εn = |ξn+1 − ξn| / |ξn+1| |
---|---|---|
".(100-$halt)." | ".$xi." | "; $epsilon = $xi; $xi = $xi - ((val($a, $xi))/(val($ap, $xi))); $epsilon = abs($epsilon - $xi)/abs($xi); echo "".$epsilon." |
Nie uzyskano zbieżności po 100 iteracjach. Metoda najprawdopodobniej nie jest zbieżna.
"; } else { $x[$j] = $xi; $xc++; } } // koniec udokładniania pierwiastków echo "