$p_0(x)=f(x)$
$p_1(x)=f'(x)$
$p_2(x)=-rem(p_0, p_1)=p_1 q_0 - p_0$
$p_3(x)=-rem(p_1, p_2)=p_2 q_1 - p_1$
...
$p_m(x)=-rem(p_{m-1}, p_{m-2}) = p_{m-1} q_{m-2} - p_{m-2}$
$rem$ oznacza resztę z dzielenia wielomianów.
Przykładowo jeśli obliczeniach otrzymamy: ++--, to $\sigma=1$, jeśli +-+- to $\sigma=2$.
Twierdzenie Sturma działa tylko dla pierwiastków pojedynczych. Jeśli wielomian posiada podwójne pierwiastki musimy się ich najpierw pozbyć budując wielomian bez nich: $\displaystyle ^P/_{GCD(P, P')}$. Przez dzielenie rozumiemy dzielenie wielomianów. GCD to największy wspólny dzielnik. W moim przypadku sekwencje tą wykorzystujemy do obliczeń związanych z raytracingiem. Tak więc póki co nie musiałem się przejmować przypadkami wielokrotnych pierwiastków. Jedynym jest torus zredukowany do sfery.
Tak więc zakładamy, że nasz wielomian posiada tylko pojedyncze pierwiastki.
Dla wielokrotnych twierdzenie nie będzie działało dobrze. Zmiana w znaku zostanie wykryta zawsze dla pierwiastka nieparzystego. Ale może zostać także stwierdzona dla parzystego. Przynajmniej tak to działa do wielomianów stopnia czwartego - sprawdzone doświadczalnie. Innymi słowy nigdy nie dostaniemy fałszywej liczby pierwiastków ale możemy nie dostać prawdziwej. W przypadku torusa zredukowanego do sfery problemem są raczej błędy numeryczne niż ten błąd.
Sekwencja Sturma posiada następujące właściwości:
$1^\circ \quad p_m(x) \neq 0$ \: dla $x \in (a, b)$
$2^\circ \quad \text{jeśli} \:\: 0<i<m-1$ \text{to} $p_{i-1}(x) = p_{i+1}(x) = 0$ \: \text{dla} \: $x \in (a, b)$
$3^\circ \quad \text{jeśli} \:\: 0<i<m$ to $p_{i-1}(x) \neq 0$ i $p_{i}(x) \neq 0$ \: \text{dla} \: $x \in (a, b)$
$4^\circ \quad$ istnieje takie $\alpha$, że dla $p_0(x)=0$ zachodzi $sign(p_0(x-\alpha)) = sign(p_0(x+\alpha))$
$1^\circ$ Wynika z definicji reszty. Musi ona być co najmniej o stopień mniejsza od dzielnej. Stopień $p_1$ też jest o stopień niższy od $p_0$. Stopień wielomianu maleje wraz z wyrazami sekwencji. Jeśli wielomian jest stopnia zerowego to nasza sekwencja liczy tylko jeden wyraz. Dla wielomianu pierwszego stopnia jego pochodna zawsze będzie różna od zera. Dla wielomianów wyższych stopni, ponieważ pierwiastki są pojedyncze to $p_0(x)$ i $p_1(x)$ nie mają wspólnych miejsc zerowych, reszta $p_2$ musi być różna od zera. I nie ma miejsc zerowych wspólnych z $p_1$ i $p_0$. W konsekwencji kolejne reszty nie mogą być zerowe, ich stopień z braku wspólnych miejsc zerowych będzie się zmniejszał zawsze o jeden do $p_{m-1}=ax-b$, $p_m=a$.
$2^\circ$ Sekwencja musi mieć co najmniej 3 wyrazy. Jeśli $p_1=0$ to z definicji wynika, że $p_2=p_1 q_0 - p_0 = -p_0$. Jeśli $p_i=0$, gdzie $i>1$ to z definicji na ogólny wyraz ciągu: $p_i(x)= p_{i-1} q_{i-2} - p_{i-2}$, wynika, że $p_{i+1}(x)= p_{i} q_{i-1} - p_{i-1}$, dla $p_i=0$ mamy: $p_{i+1}=-p_{i-1}$.
$3^\circ$ Wynika z rozważań z punktu $1^\circ$. Stopień wielomianu w sekwencji zawsze maleje od jeden do stałej i nie mają one wspólnych miejsc zerowych.
$4^\circ$ Wynika z definicji pochodnej. Pochodna wielomianu osiąga zera tylko w ekstremach wielomianu. Pomiędzy nimi, czyli w obszarze występowania miejsc zerowych wielomianu jest monotoniczna. O ile miejsca zerowe nie są wielokrotne, wtedy miejsce zerowe może być jednocześnie ekstremum. Dla pojedynczych miejsc zerowych wielomian zawsze przechodzi przez miejsce zerowe i zmienia znak.
Wszystkie te cztery właściwości pozwalają nam udowodnić poprawność twierdzenia. Udowodnimy, że funkcja $\sigma$ wyznaczająca ilość zmian znaków jest malejąca i maleje o jeden tylko przy przejściu przez miejsce zerowe wielomianu.
Sekwencja składa się z wielomianów, których miejsca zerowe są pojedyncze. Jeśli przechodzimy przez zero w wyrazie sekwencji wyraz ten zmienia znak. Z $3^\circ$ wynika, że taki wyraz otoczony jest przez wyrazy niezerowe. Ale nie jest nigdzie powiedziane, że inne wyrazy nie mogą być zerami.
Przypuśćmy, że miejsce zerowe nastąpiło w jednym lub więcej wyrazach różnych od $p_0$. Zgodnie z $4^\circ$ i $2^\circ$ mamy cztery kombinacje znaków przed miejscem zerowym: -++, --+, +--, ++- i po:
--+, -++, ++-, +--. tak więc liczba znaków w sekwencji nie ulega zmianie.
Przypuśćmy, że miejsce zerowe wystąpiło w wyrazie $p_0$ (może też wystąpić w innych wyrazach sekwencji). Zgodnie z $4^\circ$ $p_1$ nie zmienia znaku w otoczeniu miejsca zerowego. Same wielomian $p_0$ zmienia znak z uwagi na to, że posiada on tylko pojedyncze miejsca zerowe. Jeśli wielomian $p_0$ zmienia znak z + na - to znak pochodnej to -, tak więc możliwa kombinacja znaków przed to: +-, -+, zaś po: --, ++.
Widzimy więc, że tylko podczas przejścia przez miejsca zerowe wielomianu $p_0$ liczba zmian znaków w sekwencji maleje o jeden. Możemy więc powiedzieć, że liczba zmian znaków w przedziale $(a, b)$ równa się $\sigma(a) - \sigma(b)$. Ponieważ zmiana znaków następuje przy przejściu przez miejsce zerowe wielomianu $p_0$ jest to tym samym liczba miejsc zerowych w tym przedziale.
Twierdzenie Sturma wykorzystuje się do znajdowania miejsc zerowych wielomianów. Dla wielomianów stopnia 5 i wyższego nie istnieje rozwiązanie dane wzorami wprost. Poza tym wzory na pierwiastki wielomianu stopnia 3 i 4 są niestabilne numerycznie (dokładność double jest niewystarczająca by policzyć je dokładnie, lub by policzyć je dobrze kiedy są blisko siebie - wtedy często wyliczamy jako pierwiastek ekstremum między nimi). Najpierw metodą np. binarnego poszukiwania wyodrębniamy obszary w których zmiana znaku jest pojedyncza, a następnie jedną z metod aproksymacji (Newton, ModifedRegulaFasli, Bisection) przybliżamy się do miejsca zerowego, korzystając z faktu, że jest pojedyncze, i znak z jednej jego strony jest przeciwny do znaku z drugiej strony.
W przypadku obliczeń do raytracingu, interesuje nas tylko pierwsze rozwiązanie większe od zera, a ściśle od
Constants.MINIMUM_DISTANCE.Kod klasy implementującej sekwencje Sturma:
public struct SturmSequence
{
private readonly Polynomial[] m_seq;
public readonly int Length;
public SturmSequence(Polynomial a_poly)
{
a_poly = a_poly.RemoveZeroes();
m_seq = new Polynomial[a_poly.Order + 1];
Length = 1;
m_seq[0] = a_poly;
if (a_poly.Order == 0)
return;
Length++;
m_seq[1] = a_poly.Differentiate();
while (m_seq[Length - 1].Order != 0)
{
Length++;
m_seq[Length - 1] = -(m_seq[Length - 3] % m_seq[Length - 2]);
if (m_seq[Length - 1].Order == 0)
return;
}
}
public static Polynomial SingleRoots(Polynomial a_poly,
double a_precision = Constants.DOUBLE_PRECISION)
{
return a_poly / Polynomial.GCD(a_poly, a_poly.Differentiate(), a_precision);
}
public Polynomial this[int a_index]
{
get
{
if (a_index >= Length)
throw new IndexOutOfRangeException();
return m_seq[a_index];
}
}
public int NumberOfSignChangesAt(double a_x)
{
int sc = 0;
double c1 = m_seq[0].Evaluate(a_x);
for (int i = 1; i < Length; i++)
{
double c2 = m_seq[i].Evaluate(a_x);
if ((c1 * c2 < 0) || (c1 == 0.0))
sc++;
c1 = c2;
}
return sc;
}
public int NumberOfSignChangesAtPositiveInfinity()
{
int sc = 0;
double c1 = m_seq[0].LeadCoefficient;
for (int i = 1; i < Length; i++)
{
double c2 = m_seq[i].LeadCoefficient;
if ((Math.Sign(c1) * Math.Sign(c2) < 0) || ((c1 == 0.0) && (c2 != 0.0)))
sc++;
c1 = c2;
}
return sc;
}
public int NumberOfSignChangesAtNegativeInfinity()
{
int sc = 0;
double c1 = -m_seq[0].LeadCoefficient;
if (m_seq[0].Order % 2 == 0)
c1 = -c1;
for (int i = 1; i < Length; i++)
{
double c2 = -m_seq[i].LeadCoefficient;
if (m_seq[i].Order % 2 == 0)
c2 = -c2;
if ((Math.Sign(c1) * Math.Sign(c2) < 0) || ((c1 == 0.0) && (c2 != 0.0)))
sc++;
c1 = c2;
}
return sc;
}
}
Kod klasy solvera wykorzystującego powyższą sekwencje znajdujący tylko jedno miejsce zerowe większe od Constants.MINIMUM_DISTANCE większe od zera:
public static class Constants
{
public const int STURM_ANY_ROOTS_MAXIMUM_ITERATIONS = 32;
public const int STURM_FIND_ONE_ROOT_AREA_MAXIMUM_ITERATIONS = 800;
public const double MINIMAL_DISTANT = 1e-8;
}
public struct PolynomialSolverAprox
{
private readonly SturmSequence m_sturm_sequence;
private RootFinder m_root_finder;
public PolynomialSolverAprox(Polynomial a_polynomial, RootFinder a_root_finder)
{
m_sturm_sequence = new SturmSequence(a_polynomial);
m_root_finder = a_root_finder;
}
public double FindOneRootInRange(double a_a, double a_b, int a_num_roots_at_a)
{
for (int i = 0;
i < Constants.STURM_FIND_ONE_ROOT_AREA_MAXIMUM_ITERATIONS; i++)
{
double mid = (a_a + a_b) / 2;
int atmid = m_sturm_sequence.NumberOfSignChangesAt(mid);
int nroots = a_num_roots_at_a - atmid;
if (nroots == 1)
{
Polynomial p = m_sturm_sequence[0];
return m_root_finder.FindRoot(a_a, mid, (x) => p.Evaluate(x));
}
else if (nroots == 0)
a_a = mid;
else
a_b = mid;
}
return a_a;
}
public double SolveOne()
{
int nroots_at_a = m_sturm_sequence.NumberOfSignChangesAt(
Constants.MINIMAL_DISTANT);
int nroots_at_pos_inf =
m_sturm_sequence.NumberOfSignChangesAtPositiveInfinity();
int nroots = nroots_at_a - nroots_at_pos_inf;
if (nroots == 0)
return Double.PositiveInfinity;
double a = Constants.MINIMAL_DISTANT;
double b = a + 1.5;
int nroots_at_b = nroots_at_a;
for (int its = 0; its < Constants.STURM_ANY_ROOTS_MAXIMUM_ITERATIONS; its++)
{
nroots_at_b = m_sturm_sequence.NumberOfSignChangesAt(b);
if (nroots_at_a - nroots_at_b >= 1)
break;
b *= 2.0;
}
if (nroots_at_a == nroots_at_b)
return Double.PositiveInfinity;
return FindOneRootInRange(a, b, nroots_at_a);
}
}
Brak komentarzy:
Prześlij komentarz