Wissenschaftliches Rechnen

Definition: Sei \(\varepsilon \in [\tilde{\nu}, 1 – \tilde{\nu}]\). Seien nach der Trapezregel \({}^{T}{\uparrow}_{z\in A}{f(z){\downarrow}z:={\LARGE{\textbf{+}}}_{z\in A}{\tilde{2}(f(z)+f(\overset{\rightharpoonup}{z}))(\overset{\rightharpoonup}{z}-z)}}\) und nach der Mittelpunktregel \({}^{M}{\uparrow}_{z\in A}{f(z){\downarrow}z:={\LARGE{\textbf{+}}}_{z\in A}{f(\tilde{2}(z\,+\overset{\rightharpoonup}{z}}))(\overset{\rightharpoonup}{z}-z)}\). Für \(j, k \in {}^{\omega}\mathbb{N}, f \in \mathcal{C}_{\pi}^{j+2}\) und die Fourierkoeffizienten \(c_k := \tilde{\hat{\pi}}{\uparrow}_{-\pi}^{\pi}{f(t)\tilde{\epsilon}^{\underline{k}t}{{\downarrow}t}}\) heißt die Reihe \({\LARGE{\textbf{+}}}_{k=-\omega}^{\omega}{c_k\epsilon^{\underline{k}t}}\) Fourierreihe mit evtl. anderen Periodenlängen als \(\hat{\pi}\)¹. Seien \(f_n^*(z) = f(\eta_nz)\) die Schwestern zur TR \(f(z) \in \mathcal{O}(D)\) um 0 auf dem Gebiet \(D \subseteq {}^{\omega}\mathbb{C}\) mit \(m, n \in {}^{\omega}\mathbb{N}^{*}\) und \(\eta_n^m := \underline{1}^{2^{\lceil m/n \rceil}}\) sowie \(\delta_n^*f = \check{f} – \check{f}_n^*\) die halben Schwesterabstände von \(f.\) Seien \(s \in {}^{\omega}\mathbb{Z}^{*}, \mu_s^m := z^{-s}m!/(m – s)!\) und \(\widetilde{|s|\text{-}!} := 0\). Auf Ebene der TRn bilden \(\mu\) und \(\eta\) den \(\mu\)-\(\eta\)-Kalkül.\(\triangle\)

Bemerkung: Letzterer erlaubt Integrale und Ableitungen einfach und endlich geschlossen darzustellen². So gilt für die TRn \({}^2f(x){\uparrow}_0^x{\uparrow}_0^{y}{g(v){\downarrow}v{\downarrow}y} = f(x\mu_2)g(x\mu_{-2}).\)

Beispiel: Aus \(\text{Re} \; c \in [\tilde{\nu}, 1 + \tilde{\nu}[, c \in {}^{\omega}\mathbb{C}\) und \(\varepsilon := \tilde{2}^j, j \in {}^{\omega}\mathbb{N}^{*}\) folgt \(\zeta (n+c)=2^{jn}\tilde{n}{\LARGE{\textbf{+}}}_{k=1}^{n}{\delta_n^* v_c(\tilde{2}^j u^k)}\) für \(z \in {}^{1-\tilde{\nu}}\dot{\mathbb{C}} \subset \mathbb{D}\)³ und \(v_c(z):={\LARGE{\textbf{+}}}_{m=1}^{\omega }{\zeta (m+c){{z}^{m}}}=z{\LARGE{\textbf{+}}}_{m=1}^{\omega }{{{{\widetilde{m}}}^{c}}\widetilde{z-m}}.\)

Beispiel: Eine Richardson-Extrapolation bestimmt die Digammafunktion \(\psi\) und\[\zeta(\grave{n}) =\widetilde{\varepsilon^n\hat{n}}{\LARGE{\textbf{+}}}_{k=1}^n{(\psi(\varepsilon u^k \underline{1}^{2\tilde{n}}) – \psi(\varepsilon u^k))} + \mathcal{O}(\varepsilon^{\hat{n}})\] für \(n = 2\) und \(\varepsilon = \tilde{2}^{13}\) als \(\zeta(3) = 2^{24}{\LARGE{\textbf{+}}}_{k=1}^2{(\psi(\underline{\varepsilon} u^k) – \psi(\varepsilon u^k))} + \mathcal{O}(10^{-16})\).

Darstellungssatz für Integrale: Die TR (s. u.) \(f(z) \in \mathcal{O}(D)\) um 0 für \(D \subseteq {}^{\omega}\mathbb{C}\) ergibt mit \(\grave{m}, n \in {}^{\omega}\mathbb{N}^*\)\[{\uparrow}_0^z…{\uparrow}_0^{\zeta_2}{f(\zeta_1){\downarrow}\zeta_1\;…\;{\downarrow}\zeta_n} = f(z\mu_n).\square\]Darstellungssatz für Ableitungen: Mit \({}^{\widetilde{\nu}}\dot{\mathbb{C}} \subset \mathbb{D} \subseteq {}^{\omega}\mathbb{C}\) ergeben die \(n\)-ten Einheitswurzeln, die TR\[f(z):=f(0) + {\LARGE{\textbf{+}}}_{m=1}^{\omega }{\widetilde{m!}\ {}^mf(0){z^m}},\]\(\varepsilon := \tilde{2}^j\tilde{r}, j \in {}^{\omega}\mathbb{Z}, n = \epsilon^{\sigma} \in {}^{\omega}\mathbb{N}^{*}, u :=\epsilon^{\tilde{n} \hat{\underline{\pi}}}\) und der Konvergenzradius \(r \in {}^{\nu}{\mathbb{R}}_{>0}\) von \(f\)\[{}^nf(0)=2^{jn}\acute{n}!{\LARGE{\textbf{+}}}_{k=1}^{n}{\delta_n^* f(\tilde{2}^j u^k)}.\square\]Universeller Mehrschrittsatz: Mit \(n \in {}^{\nu}\mathbb{N}_{\le p}, k, m, p \in {}^{\nu}\mathbb{N}^{*}, {\downarrow}\overset{\rightharpoonup}{x} \in\, ]0, 1[, x \in [a, b] \subseteq {}^{\omega}\mathbb{R}, y : [a, b] \rightarrow {}^{\omega}\mathbb{R}^q,\) \(f : [a, b] \times {}^{\omega}\mathbb{R}^{q \times n} \rightarrow {}^{\omega}\mathbb{R}^q, g_k(\overset{\rightharpoonup}{x}) := g_{\acute{k}}(x)\) und \(g_0(a) = f(\overset{\leftharpoonup}{a}, y_0, … , y_{\acute{n}})\) ergibt die TR des Anfangswertproblems \(n\)-ter Ordnung \(y^\prime(x) = f(x, y((\rightharpoonup)^0 x), … , y((\rightharpoonup)^{\acute{n}} x))\)\[y(\overset{\rightharpoonup}{x}) = y(x) + {\downarrow}\overset{\rightharpoonup}{x}{\LARGE{\textbf{$\pm$}}}_{k=1}^{p}{\left (g_{p-k}(\overset{\rightharpoonup}{x}){\LARGE{\textbf{+}}}_{m=k}^{p}{\widetilde{m!}\tbinom{\acute{m}}{\acute{k}}}\right )} + \mathcal{O}(({\downarrow}\overset{\rightharpoonup}{x})^{\grave{p}}).\square\]Auf- und Ableitungssatz für eine KR: Mit \(m \in {}^{\omega}\mathbb{Z}\), \(q = \tilde{\varepsilon}(z-a)\), \(a \in \mathbb{D}\) und \(j, k \in \mathbb{N}_{<n}\) ergeben modulare Arithmetik⁴ und die \(n\)-ten Einheitswurzeln zur analogen TR:\[{}^mf_{\tilde{n}}(z) := \tilde{n}(q^j)^T(\delta_{jk}\widetilde{\varepsilon q}^mj!/(j-m)!)({\tilde{u}}^{jk})(f({\varepsilon u}^k+a))+R_n^{m}(z).\square\]Umkehrsatz für KRn: Für \(y \in f(\mathbb{D}), y(a) = b\) und \({}^1y(a) \ne 0\) ergibt der Satz von Bürmann⁵:\[f_{\tilde{n}}^{-1}(y) := a + \tilde{n} {\LARGE{\textbf{+}}}_{m=1}^n{\widetilde{m}{\tilde{\varepsilon}}^{\acute{m}}(y – b)^m({\tilde{u}}^{\acute{m}k})^T(f(\varepsilon u^k + a)^{-m})}+R_n^{-}(y).\square\]Bemerkung: Ein hinreichend kleines \(\tilde{n}|a – z_0|\) liefert für \(y = 0\) eine Nullstelle \(z_0 \in {}^{\omega}\mathbb{C}\) von \(f(z) \in {}^{\omega}\mathbb{C}\). Der vorige Satz ermöglicht implizite Differenzialgleichungen für KRn einfach in explizite umzuwandeln.

Ableitungssatz für Fourierreihen: Mit \(f \in \mathcal{C}_{\pi}^{m+2}\)⁶, \(m \in {}^{\omega}\mathbb{N}, j \in [-\acute{n}, n] \cap \mathbb{Z}, t \in [-\pi, \pi]\) und \(k \in \mathbb{N}_{<\hat{n}}\) ergibt sich mit der Möglichkeit zur gliedweisen Integration folgendes KR-Analogon:\[{}^m\check{f}_{\tilde{n}}(t):=\tilde{n}(u^{\tilde{\pi}jnt})^T(\delta_{(j+\acute{n})k}\underline{j}^m)(\tilde{u}^{jk})(\check{f}(\tilde{n}\pi k – \pi))+\mathcal{O}(\tilde{n}).\square\]Folgerung: Die Stützstellen \(kr := \tilde{n}\pi k\) des glatten \(f(kr)\) liefern mit \(j\) wie \(k\) interpolierend in \(\mathcal{O}(\sigma n)\):\[{}^m\check{f}_{\tilde{n}}(t):=\tilde{n}(u^{\tilde{r}jt})^T(\delta_{jk}\underline{j}^m)(\tilde{u}^{jk})(\check{f}(kr)).\square\]Bemerkung: Die Identität anstatt \(\delta_n^*\) liefert beliebig genaue Näherungen für die \({}^nf\). Der Fehler für analog definierte \(m\)-dimensionale KRn mit je \(\binom{m+n}{n}\) Ableitungen ist vergleichbar.

Satz (KRn-Fixpunktverfahren): Jedes \(z \in {}^{\omega}\mathbb{C}\) eines beliebigen \(m\)-Polynoms \(p(z) = 0\) mit \(m \in [2, n] \cap \mathbb{N}\) für \(n := 2^r, r \in {}^{\nu}\mathbb{N}^*\) und Koeffizienten aus \({}^{\nu}\mathbb{C}\) lässt sich (als Anschub) ebenfalls in \(\mathcal{O}(\sigma n)\) bestimmen.

Beweis und Algorithmus: Sei \(U = (\tilde{u}^{jk})\) mit \(j, k \in \mathbb{N}_{<n}, u :=\epsilon^{\tilde{n} \hat{\underline{\pi}}}, q := \hat{z}\) und \(s_k := p(\tilde{2}u^k)\). Eine einfache Transformation erreicht \(|q| < \tilde{2}\) für alle Nullstellen \(\zeta\) von \(p(z)\) und \(p(0) = 1\). Für \(p(z) = \tilde{n}(q^j)^TUs = \tilde{n}\mu^Ts = 0\) ergibt sich die vereinfachte Iteration \(\mu^* = U_{1}^{-T}\mu U((\delta_{jk}\tilde{u}^j)U^{-1}\mu-(U_{\acute{n}}^{-T}\mu+\beta s^T\mu, \beta s^T\mu, …, \beta s^T\mu)^T)\). Hierbei sei \(\delta_{jk}\) das Kronecker-Delta, der Startpunkt \(q := \tilde{2}\) und \(\beta \in {}^{\nu}\mathbb{C}^*\), sodass jeweils \(||\mu^*-\mu||\) sich ungefähr halbiert und \(\mu^Ts = 0\) gilt. Für \(m > 2\) bildet Polynomdivision den Schluss.\(\square\)

Approximationssatz: Für \(x \in {}^{\omega}\mathbb{R}\) und \(m, n \in {}^{\nu}\mathbb{N}\) erlauben die \({}^sf(x) \in {}^{\omega}\mathbb{R}\) wie oben die Interpolante

\(g(x) := {\LARGE{\textbf{+}}}_{r=0}^{\acute{m}}{\chi_{]x_r, x_{\grave{r}}[}(x)((x_{\grave{r}}-x)p_r(x)+(x-x_r)p_{\grave{r}}(x))/(x_{\grave{r}}-x_r)}+\)\({\LARGE{\textbf{+}}}_{r=0}^m{\chi_{\{x_r\}}(x)p_r(x)}\) mit \(p_r(x) := {\LARGE{\textbf{+}}}_{s=0}^n{{}^sf(x_r){(x-x_r)}^s/s!}\)

in \(\mathcal{O}(\sigma mn)\) zu berechnen, wobei \({}^sf(x_r) = {}^sg(x_r)\) für alle \(x_r \in {}^{\omega}\mathbb{R}\) gilt. Im Komplexen ist \({}^{\omega}\mathbb{R}\) durch \({}^{\omega}\mathbb{C}\) zu ersetzen und es gelte \(x = \gamma(t) \in {}^{\omega}\mathbb{C}\) für den Weg \(\gamma(t)\) mit \(t \in {}^{\omega}\mathbb{R}.\square\)

Differenzialgleichungssatz: Sei \(d(z) \in {}^{\omega}\mathbb{C}^{n+2}\) Vektor von KRn mit \(n \in {}^{\omega}\mathbb{N}\). Die Lösung \(y_0\) des Systems \({}^{\grave{n}}y = \left(\complement_0^n\ {}^{m}y, 1\right)d(z) \in {}^{\omega}\mathbb{C}\) mit / ohne \(\left(\complement_0^n\ {}^{m}y(a)\right)^T\in {}^{\omega}\mathbb{C}^{\grave{n}}\) und \(a, z \in {}^{\omega}\mathbb{C}\) lässt sich ebenso in \(\mathcal{O}(n^3)\) bestimmen.\(\square\)

Bemerkung: Werden Funktionswerte von Kronpunkten \(\varepsilon u^m + a\) um den Kronmittelpunkt \(a\) ausgewertet, wobei \(\varepsilon = \tilde{2}\) und \(n = 64\) gut gewählt sind, und abbrechende KRn (evtl. auch Laurent- und Fourierreihen) genutzt, lassen sich sowohl eine Polynomdivision ausführen als auch (zumindest abschnittsweise) analytische, nichtlineare bzw. partielle (Integro-)Differentialgleichungssysteme lösen.

Beispiel: Mit \(b = (f(\varepsilon u^0 + a), …, f(\varepsilon u^{\acute{n}} + a), 1)^T \in {}^{\omega}\mathbb{C}^{\grave{n}}\) und einer oberen Dreiecksmatrix \(T(z) \in {}^{\omega}\mathbb{C}^{\grave{n}\times\grave{n}}\) haben gewöhnliche Differenzialgleichungen nur Lösungen, wenn \(T(z)b = 0\) gilt.

Reihensatz für Integrale: Für \(c \in [0, a], f \in {}^{\nu}\mathcal{C}^{n}\) und \(n, \hat{p}, r\,(= 2^p), t\,(=2^v), \hat{v} \in {}^{\nu}2\mathbb{N}^*\) zeigen KRn

\({\uparrow}_0^a f(w){\downarrow}w = {\uparrow}_0^{\tilde{t}a}{\LARGE{\textbf{+}}}_{s=1}^t{f(x+\acute{s}\tilde{t}a){\downarrow}x} = {\uparrow}_0^1 g(y){\downarrow}y =\)\( {\LARGE{\textbf{+}}}_{\check{q}=1}^{\check{r}}{{\LARGE{\textbf{+}}}_{\check{m}=0}^{\check{n}-1}{\widetilde{\grave{m}!}\ {}^m\hat{g}\left (\tilde{r}\acute{q} \right){\tilde{r}}^{\grave{m}}}}+\mathcal{O}\left(\widetilde{\grave{n}!}\ {}^ng(\tilde{a}c){\tilde{r}}^n\right).\square\)

Bemerkungen: KRn-Äquivalente können die \({}^mg\) mit \(g(y) := af(ay)\) und \(a > 0\) ersetzen. Die Transformation \(z := \epsilon^{\pm x}\) macht hier unendliche Integralgrenzen endlich. Ein endliches Integral verringert den Betrag des Restglieds hinreichend. Die Mittelpunktregel ist hier vorteilhafter als die Trapezregel.

Beispiel: Mit \({}_{\mp}g(\vartheta) := \check{\pi}(1 – {\varepsilon}^2 {\sin}^2(\check{\pi}\vartheta))^{\mp\tilde{2}}\) ist das vollständige elliptische Integral I. und II. Art\[{\uparrow}_0^{1} {}_{\mp}g(\vartheta) {\downarrow}\vartheta = {}_{\mp}g(\check{1}) + \widetilde{24}\ {}_{\mp}^2g(\check{1}) + \widetilde{1920}\ {}_{\mp}^4 g(\check{1}) + \mathcal{O}\left(\widetilde{9!}\ {}_{\mp}^6g(\check{1})\right).\]Zweite Euler-Maclaurin-Formel: Für \(f(\check{q}) = g(\tilde{r}\acute{q})\) und \(k = \grave{a} = \grave{r}\) liefert der vorangehende Satz in \(\mathcal{O}(\sigma n)\)\[{\LARGE{\textbf{+}}}_{\check{q}=1}^{\check{r}} f(\check{q}) = {\uparrow}_{\check{1}}^{\check{k}} f(x){\downarrow}x + {\LARGE{\textbf{+}}}_{\check{m}=1}^{\check{n}-1} H_m \left ({}^{\acute{m}}f(\check{k}) – {}^{\acute{m}}f(\check{1}) \right ) + \mathcal{O}\left (H_n \left ({}^{\acute{n}}f(\check{k}) – {}^{\acute{n}}f(\check{1}) \right )\right ).\]Beweis: Für \(h(x) = x/\)sin \(x\) und \(H_m := \underline{1}^m\widetilde{m!}\ {}^mh(0)=\widetilde{m!}B_m(2-2^m) \rightarrow 2\tilde{\underline{\pi}}\text{-}^m\) folgt aus\[{\LARGE{\textbf{+}}}_{\check{q}=1}^{\check{r}} g(\tilde{r}\acute{q}) = \check{r}{\uparrow}_0^1 g(y){\downarrow}y – {\LARGE{\textbf{+}}}_{\check{q}=1}^{\check{r}}{\LARGE{\textbf{+}}}_{\check{m}=1}^{\check{n}-1} {\widetilde{\grave{m}!}\ {}^mg(\tilde{r}\acute{q})} + \mathcal{O}\left(\widetilde{\grave{n}!}\ {}^ng(\tilde{a}c){\tilde{r}}^n\right)\]durch Einsetzen und Zusammenfassen der in Abhängigkeit von Partitionen gezählten Fakultäten die Behauptung, wobei die \(B_m\) (vgl. Nichtstandardanalysis) Bernoulli-Zahlen sind und \(B_{\grave{m}} = 0, B_0=1\) sowie \(B_1=-\tilde{2}\) gilt.\(\square\)

Erste Euler-Maclaurin-Formel: Vollständige Induktion nach \(n\) zeigt ebenso⁷

\({\LARGE{\textbf{+}}}_{\check{q}=0}^{\check{r}}f(\check{q})=\uparrow_0^{\check{r}}{f\left(x\right){\downarrow}x}+\check{f}(\check{r})+\check{f}(0)+{\LARGE{\textbf{+}}}_{\check{m}=1}^{\check{n}-1}{\widetilde{m!}}B_m\left({}^{\acute{m}}f(\check{r})-{}^{\acute{m}}f(0)\right)+\)\(\mathcal{O}\left({\widetilde{n!}}B_n\left({}^{\acute{n}}f(\check{r})-{}^{\acute{n}}f(0)\right)\right).\square\)

Bemerkung: Ersetze die Ableitungen der Euler-Maclaurin-Formeln auch durch KRn-Äquivalente.

Das Kron-FFT-Aberth-Verfahren zur Nullstellenbestimmung

Sei

Gesucht sind die Nullstellen von \(p\). Im Gegensatz zu einem lokalen Newton-Verfahren wird nicht vorausgesetzt, dass bereits ein Startwert in der Nähe einer einzelnen Nullstelle bekannt ist. Stattdessen wird das Polynom zuerst durch eine Kronverschiebung global um einen frei wählbaren Mittelpunkt dargestellt.

Sei \( a\in{}^{\nu}\mathbb C \) ein Kronmittelpunkt. Weiter seien \( n=2^r, r\in{}^{\nu}\mathbb N^*, n>m, \) und \( u:=\epsilon^{\tilde n\hat{\underline{\pi}}} \) eine primitive \(n\)-te Einheitswurzel. Für \(\sigma\in{}^{\nu}\mathbb Z\) sei der Kronradius \(\rho:=2^\sigma\). Diese Wahl ist numerisch günstig, da die Skalierung mit Potenzen von \(\rho\) binär exakt erfolgt, solange kein Über- oder Unterlauf auftritt.

Die Kronpunkte lauten

Setze

Skalierte KR

Mit der skalierten Variable \( z=a+\rho q \) gilt

Die Koeffizienten \(d_j\) ergeben sich aus den Kronwerten durch Einheitswurzelprojektion:

Die unskalierte TR beziehungsweise KR in \( h=z-a \) mit \( b_j=\tilde\rho^{\,j}d_j \) lautet

Für die Nullstellensuche ist die skalierte Form \(P(q)\) günstiger, da die gesuchten \(q\)-Nullstellen bei geeigneter Wahl von \(\rho\) in einem moderaten Bereich, häufig nahe dem Einheitskreis, liegen.

FFT-Fassung

Die Folge der Kronwerte \( s=(s_0,\dots,s_{\acute n})^\top \) wird per schneller Fourier-Transformation projiziert. Bis auf die Vorzeichenkonvention der verwendeten FFT gilt

Für Polynome mit \(m<n\) ist die Rekonstruktion im exakten Körper aliasfrei. Numerisch empfiehlt sich \( n\ge \widehat{m}+1, \) meist wiederum als Zweierpotenz.

Initialwerte für Aberth-Ehrlich

Für das skalierte Polynom

\[ P(q)= {\LARGE{\textbf{+}}}_{j=0}^{m}{d_jq^j} \]

werden \(m\) verschiedene Anfangswerte benötigt. Sie werden nicht aus der Nähe einzelner Nullstellen gewonnen, sondern gleichmäßig auf einem Startkreis gewählt.

Sei \( v:=\epsilon^{\tilde m\hat{\underline{\pi}}} \) eine primitive \(m\)-te Einheitswurzel und sei \( R_A\in{}^{\nu}\mathbb R_{>0} \) ein Aberth-Startkreisradius. Dann werden die Anfangswerte gesetzt als

In gewöhnlicher Schreibweise ist dies

Der halbe Winkelschritt \(i-\check 1\) vermeidet einfache Symmetrie-Sonderfälle. Wird \(\rho\) so gewählt, dass die relevanten Nullstellen von \(p\) näherungsweise durch \( |z_\ell-a|\le \rho \) erfasst werden, dann liegen die zugehörigen skalierten Nullstellen \( q_\ell=\tilde\rho(z_\ell-a) \) typischerweise im Einheitskreis. Dann ist \( R_A=1 \) eine natürliche Wahl.

Ist keine geometrische Information verfügbar, lässt sich eine Cauchy-Schranke für \(P\) verwenden:

\[ \acute{R}_A = |\tilde{d}_m|\max_{0\le j<m}|d_j|, \qquad d_m\ne0. \]

Dann liegen alle Nullstellen von \(P\) in der Kreisscheibe \( |q|\le R_A. \)

Aberth-Ehrlich-Schritt ohne Eigenwertzerlegung

Für aktuelle Näherungen \( q_1,\dots,q_m \) sei

Dann lautet der Aberth-Ehrlich-Schritt⁸

Der Newton-Anteil \(\Delta_i\) zieht \(q_i\) zu einer Nullstelle von \(P\). Der zweite Term im Nenner wirkt als gegenseitige Abstoßung der Näherungen. Dadurch wird verhindert, dass mehrere Näherungen ohne Weiteres auf dieselbe Nullstelle kollabieren. Nach der Iteration erfolgt die Rückskalierung \(z_i=a+\rho q_i.\)

Abbruchkriterium

Ein rein schrittweitenbasiertes Abbruchkriterium ist bei eng beieinanderliegenden Nullstellen oder Clustern nicht ausreichend. Daher wird zusätzlich ein skalierter Auswertungsfehler verwendet:

Nach einer Mindestzahl von Iterationen wird abgebrochen, wenn zugleich

und

gelten. Dadurch wird ein vorschneller Abbruch bei numerisch empfindlichen Clustern vermieden.

Politur im Anwendungsmodus

In einer realen Anwendung ist \(p\) entweder durch Koeffizienten, durch eine Auswertungsvorschrift oder durch eine lokal rekonstruierte KR gegeben. Mit \(p_1(z_i) := {}^1p(z_i) \ne 0\) darf die abschließende Politur nur an dieser tatsächlich verfügbaren Darstellung erfolgen:

Bei Koeffizienten erfolgt die Auswertung beispielsweise mit dem Horner-Schema. Bei einer Black-Box-Funktion wird die Originalfunktion ausgewertet.

Mehrfachnullstellen

Mehrfachnullstellen sind intrinsisch schlecht konditioniert. Eine kleine Störung der Koeffizienten kann eine \(r\)-fache Nullstelle in \(r\) nahe einfache Nullstellen aufspalten. Dies ist keine Schwäche der Kronverschiebung, sondern eine Eigenschaft der Nullstellenabbildung.

Sind Koeffizienten mit hinreichender Genauigkeit verfügbar, ist eine quadratfreie Zerlegung zweckmäßig. Sei \( g:=\gcd(p,{}^1p). \) Dann ist \( p_{\mathrm{sf}}:=p/g \) der quadratfreie Anteil von \(p\). Die Nullstellen von \(p_{\mathrm{sf}}\) sind einfach. Die Vielfachheiten werden anschließend aus den abgespaltenen Faktoren rekonstruiert.

Ist eine Vielfachheit \(r\) lokal bekannt oder numerisch geschätzt, so lässt sich statt des einfachen Newton-Schritts die modifizierte Form verwenden:

Präzisionsadaptive Fassung

Bei Clustern, nahezu mehrfachen Nullstellen oder schlecht skalierten Polynomen werde die Arbeitsgenauigkeit erhöht. Sei \(\beta\) die Zahl der verwendeten Bits. Wird das Residualkriterium nicht erreicht oder bleibt ein Cluster instabil, so wird \(\beta\leftarrow\hat\beta\) gesetzt, also die Arbeitsgenauigkeit verdoppelt. Die Kronauswertung, die Fourier-Projektion, der Aberth-Ehrlich-Schritt und die Politur werden dann mit erhöhter Genauigkeit wiederholt. Die Skalierung \( z=a+\rho q \) mit \( \rho=2^\sigma \) bleibt dabei erhalten.

Benchmarkmodus

Für synthetische Tests lassen sich Nullstellen \(r_1,\dots,r_m\) vorgeben und das Testpolynom \[ p(z)={\LARGE{\textbf{$\times$}}}_{j=1}^{m}{(z-r_j)} \] faktorisiert auswerten. Diese faktorisierte Auswertung vermeidet die schlecht konditionierte Rückbildung von Nullstellen auf Monomkoeffizienten. Insbesondere wird dadurch der Wilkinson-Effekt der Testgenerierung nicht mit dem Fehler des Nullstellenverfahrens verwechselt.

Eine im Benchmarkmodus mögliche Politur an den bekannten Faktoren \[ p(z)={\LARGE{\textbf{$\times$}}}_{j=1}^{m}{(z-r_j)} \] ist nur ein Testorakel. Sie dient dazu, die Güte der von Kron-FFT und Aberth-Ehrlich erzeugten Vorabschätzungen zu isolieren. Sie ist nicht Bestandteil des realen Lösers, wenn die Nullstellen \(r_j\) unbekannt sind.

Algorithmus

Begründung

Aus

folgt mit der Orthogonalität der Einheitswurzeln

Daher gilt

Mit \(z=a+\rho q\) ist die skalierte KR \(P(q)=p(z)\) exakt rekonstruiert, solange \(m<n\) gilt. Da dies eine bijektive affine Transformation ist, gilt

Die Aberth-Ehrlich-Iteration liefert simultane Näherungen für die Nullstellen von \(P\). Die Rückskalierung liefert daraus Näherungen für die Nullstellen von \(p\). Die Politur am ursprünglichen \(p\) vermindert nur den numerischen Restfehler und ändert das Prinzip nicht.

Aufwand

Sei \(C_p\) der Aufwand einer Auswertung von \(p\). Die Kronauswertung kostet \(\mathcal O(nC_p),\) die FFT \(\mathcal O(n\;{}_2n).\) Für \(\ell\) Aberth-Ehrlich-Iterationen ergibt sich \(\mathcal O(\ell m^2)\) bei Grad \(m\), da die gegenseitigen Abstoßungsterme quadratisch viele Differenzen enthalten. Für \(n=\mathcal O(m)\) folgt insgesamt

Bei koeffizientengegebenen Polynomen ist typischerweise \( C_p=\mathcal O(m). \) Für einfache, gut getrennte Nullstellen ist \(\ell\) in der Praxis meist moderat. Bei Clustern beziehungsweise nahezu oder tatsächlich mehrfachen Nullstellen kann \(\ell\) deutlich größer werden; dann sind angepasste Präzision oder quadratfreie Zerlegung vorzuziehen.

Klassischer Vergleich

Beim klassischen Aberth-Verfahren auf dem ursprünglichen Polynom

werden die Startwerte häufig direkt auf einem Cauchy-Kreis mit \[ \acute R= |\tilde{c}_m| \max_{0\le j<m}|c_j| \] gewählt:

Das Kron-FFT-Aberth-Verfahren verschiebt und skaliert das Problem dagegen zuerst durch \(z=a+\rho q.\) Dadurch wird die Initialisierung in der \(q\)-Ebene oft einfacher, typischerweise mit \(R_A=1\). Die Startwerte müssen nicht nahe an einzelnen Nullstellen liegen; sie müssen lediglich das relevante Nullstellengebiet angemessen umfassen.

Einordnung

Die Bestandteile des Verfahrens sind getrennt zu betrachten. Die Kron-FFT erzeugt eine verschobene globale Polynomdarstellung. Aberth-Ehrlich ersetzt die Eigenwertzerlegung der Companion-Matrix. Newton beziehungsweise modifizierter Newton poliert am ursprünglichen Problem. Das Verfahren ist daher kein lokales Newton-Verfahren mit geratenem Startwert, sondern ein globales Initialisierungsverfahren durch Kronverschiebung. Es vermeidet die Voraussetzung, bereits in der Nähe einer einzelnen Nullstelle zu starten.\(\square\)

Algorithmische Neufundierung: Abspaltung von Linearfaktoren statt blinder Auswertung

Standardmäßige Approximationsalgorithmen operieren oft strukturblind. Sie werten die Zielfunktion an diskreten Stützstellen aus oder berechnen lokale Ableitungen, ohne deren globale Topologie zu hinterfragen. Trifft ein klassischer Algorithmus auf eine Polstelle, interpretiert er diese numerisch schlicht als extrem steilen Anstieg. Das Resultat ist ein prozessualer Zusammenbruch: Lineare Gleichungssysteme werden singulär, Rundungsfehler eskalieren unkontrolliert zu NaN-Werten (Not a Number), und das resultierende Polynom verliert seine mathematische Aussagekraft.

Der hier verfolgte Ansatz verlässt diesen rein reaktiven Pfad und implementiert eine vorgezogene strukturelle Abspaltung (in Linearfaktoren). Der Berechnungsprozess gliedert sich dabei in drei strikt getrennte Phasen:

1. Struktur-Scan: Vor jeglicher Reihenentwicklung sucht ein vorgeschalteter Algorithmus aktiv nach echten Singularitäten und Nullstellen. Durch die rigorose Auswertung des Betrags komplexer Zahlen werden dabei auch Singularitäten auf der imaginären Achse fehlerfrei identifiziert.
2. Analytische Abspaltung von Linearfaktoren: Die erkannten algebraischen Strukturen werden direkt aus der Ursprungsfunktion extrahiert. Polstellen werden herausmultipliziert, Nullstellen dividiert. Übrig bleibt eine mathematisch geglättete, analytische Restfunktion, die eine numerische Eskalation vermeidet.
3. Approximation und Rekonstruktion: Erst auf diese stabile Restfunktion werden die eigentlichen Approximationswerkzeuge angewendet. Da das Verhalten der Restfunktion meist gutartig ist, genügen für die Näherung deutlich geringere Polynomgrade. Im abschließenden Rekonstruktionsschritt werden die zuvor isolierten algebraischen Pol- und Nullstellen exakt und verlustfrei in das Zähler- beziehungsweise Nennerpolynom des finalen Modells wieder integriert.

Der blinde Fleck der reellen Achse: Komplexe Nah-Singularitäten

Ein eindimensionaler Scan entlang der reellen Achse ist hochgradig effizient, besitzt jedoch einen fundamentalen blinden Fleck: komplexe Nah-Singularitäten. Es werde die rationale Funktion \(f(z) = \frac{z^2 + \epsilon^2}{z^2 + \delta^2}\) für sehr kleine, reelle Werte von \(\epsilon\) und \(\delta\) betrachtet.

Auf der reellen Achse ist diese Funktion stetig, glatt, strikt positiv und vollständig regulär. Ein rein reeller Struktursuchalgorithmus findet keine Anomalien. Dennoch verbergen sich in unmittelbarer Nähe zur reellen Achse rein komplexe Nullstellen bei \(\pm i\epsilon\) und Polstellen bei \(\pm i\delta\).

Diese verborgenen Pole fungieren als analytische Barrieren. Sie begrenzen den Konvergenzradius klassischer Potenzreihen auf \(R = \delta\), was jede globale Approximation jenseits dieses engen Intervalls zum Scheitern verurteilt. Um die Abspaltung von Linearfaktoren wahrhaft robust zu machen, darf die strukturelle Analyse daher nicht auf die mathematische Ideallinie der reellen Achse beschränkt bleiben. Sie muss zwingend zu einem zweidimensionalen Toleranzschlauch (Strip) ausgeweitet werden, um auch jene komplexen Singularitäten zu erfassen, deren Einflussgebiet bis tief in den reellen Definitionsbereich hineinragt.

Der Paradigmenwechsel – die methodische Trennung von definierender algebraischer Struktur und analytischem Rest – schützt das System auf fundamentaler Ebene vor numerischer Überanpassung (Overfitting). Die Berechnung vermeidet Singularitäten durch fehleranfällige Oszillationen zu simulieren; stattdessen passt sich das Modell von Grund auf an die wahre, innere Struktur der Funktion an.

ODE-Beispiel: Kronreihe und Faltung (FFT) vs. Runge–Kutta

Als bewusst schlichtes, aber nichtlineares Referenzproblem diene \[ {^1}y(x) = y(x)^2,\qquad y(0)=1. \] Die exakte Lösung ist \[ y(x)=\widetilde{\acute{x}\text{-}}, \] mit einer Polstelle bei \(x=1\). Damit ist der Konvergenzradius der TR um \(x_0\) gleich \(R(x_0)=\acute{x}_0\)-.

Kronreihe als Taylorentwicklung

Für einen Schritt von \(x_0\) nach \(x_0+h\) wird die lokale Reihenentwicklung \[ y(x_0+s) = {\LARGE{\textbf{+}}}_{m=0}^{n} a_m\, s^m \;+\; \mathcal{O}(s^{\grave{n}}) \] verwendet (hier ist \(s=x-x_0\)). Dann gilt \[ {^1}y(x_0+s) = {\LARGE{\textbf{+}}}_{m=0}^{\acute{n}} \grave{m}\,a_{\grave{m}}\, s^m \;+\; \mathcal{O}(s^{n}). \] Für das Quadrat entsteht das Cauchy-Produkt (Faltung) \[ y(x_0+s)^2 = \left({\LARGE{\textbf{+}}}_{m=0}^{n} a_m s^m\right)\left({\LARGE{\textbf{+}}}_{m=0}^{n} a_m s^m\right) = {\LARGE{\textbf{+}}}_{m=0}^{\hat{n}} c_m s^m \] mit \[ c_m := {\LARGE{\textbf{+}}}_{j=0}^{m} a_j\, a_{m-j}. \] Koeffizientenvergleich in \({^1}y=y^2\) liefert für \(m=0,\dots,\acute{n}\) die Rekursion \[ \grave{m}\;a_{\grave{m}} = c_m = {\LARGE{\textbf{+}}}_{j=0}^{m} a_j\, a_{m-j}, \qquad a_0 = y(x_0). \] Damit sind die \(a_m\) vollständig durch Faltungen bestimmt und enthalten keine Binomialkoeffizienten.

Faltung per FFT (Kernidee)

Für die diskrete Faltung werden Vektoren eingesetzt. Sei \(A=(a_0,\dots,a_n)\). Für eine schnelle Faltung werden Nullpolsterung auf Länge \(\acute{m} \ge \hat{n}\) (typisch \(m=2^{\lceil {}_2(\hat{n}+1)\rceil}\)) und die DFT⁹ verwendet: \[ \breve{A}_k := {\LARGE{\textbf{+}}}_{j=0}^{\acute{m}} A_j\, u^{jk}, \qquad u := \exp(\underline{\hat{\pi}}\tilde{m}), \qquad k=0,\dots,\acute{m}. \] Dann gilt für die Faltungskoeffizienten \(C=A*A\) komponentenweise \[ \breve{C}_k = \breve{A}_k^2, \qquad C_j = \tilde{m}\,{\LARGE{\textbf{+}}}_{k=0}^{\acute{m}} \breve{C}_k\, \tilde{u}^{jk}. \] Die FFT setzt diese Transformationen in \(\mathcal{O}(m\;{}_2m)\) um. Wichtig: Genau hier erweist sich die KR als praktisch, da sich Ableitungen und Produkte vektoriell — und damit FFT-tauglich — behandeln lassen.

Auswertung des Schritts

Ist \(A=(a_0,\dots,a_n)\) bestimmt, folgt der Schrittwert aus der Polynomauswertung \[ y(x_0+h)\approx {\LARGE{\textbf{+}}}_{k=0}^{n} a_k\, h^k, \] z. B. per Horner-Schema (nur Multiplikationen/Additionen).

Konkrete Rechnung für \(x_0=0\), \(h=0{,}1\), \(n=16\)

Hier ist die exakte Lösung \(y(x)=\widetilde{\acute{x}\text{-}}\), also \[ a_k = 1\quad (k\ge 0), \] denn \[ \widetilde{\acute{s}\text{-}} = {\LARGE{\textbf{+}}}_{k=0}^{\omega} s^k. \] Damit ergibt sich die \(n\)-te Partialsumme \[ y(0{,}1)\approx {\LARGE{\textbf{+}}}_{k=0}^{16} (0{,}1)^k = \frac{1-(0{,}1)^{17}}{1-0{,}1}. \] Der exakte Wert ist \[ y(0{,}1)=\frac{1}{0{,}9}=1{,}\overline{1} \] und der Restterm beträgt \[ \left|y(0{,}1)-{\LARGE{\textbf{+}}}_{k=0}^{16} (0{,}1)^k\right| = \frac{(0{,}1)^{17}}{0{,}9} \approx 1{,}111111111111111111\cdot 10^{-17}. \]

Gegenüberstellung: klassisches Runge–Kutta 4. Ordnung

Für \({^1}y=y^2\) lautet ein RK4-Schritt¹⁰ mit \(y_0=y(x_0)\), Schrittweite \(h\): \[ \begin{array}{rcl} u_1 &=& y_0^2,\\ u_2 &=& \left(y_0+h\check{u}_1\right)^2,\\ u_3 &=& \left(y_0+h\check{u}_2\right)^2,\\ u_4 &=& \left(y_0+h u_3\right)^2,\\[0.4em] y_1 &=& y_0 + \tilde{6}h\,(u_1+\hat{u}_2+\hat{u}_3+u_4). \end{array} \] Mit \(x_0=0\), \(y_0=1\), \(h=0{,}1\) ergibt sich numerisch \[ y_{\mathrm{RK4}}(0{,}1)\approx 1{,}1111104900521944, \] und damit \[ \left|y(0{,}1)-y_{\mathrm{RK4}}(0{,}1)\right| \approx 6{,}21\cdot 10^{-7}. \]

Fairer Vergleich (Prinzip)

• Genauigkeitsziel: Für hohe Genauigkeiten (z. B. \(10^{-16}\), \(10^{-34}\)) muss RK4 die Schrittweite stark verkleinern oder die Ordnung erhöhen; andernfalls dominiert der lokale Restterm.
• Kron/Faltung: Ein Schritt lässt sich (bei analytischer Lösung) deutlich größer wählen, solange \(h<R(x_0)\) bleibt, wobei primär \(n\) (Reihentiefe) und Rechenpräzision die Genauigkeit kontrollieren.
• Laufzeit: RK4 kostet pro Schritt wenige Funktionsauswertungen. Kron/Faltung erfordert pro Schritt, viele Koeffizienten (Produkte/Faltungen) plus eine Polynomauswertung zu bestimmen. Ersetzt ein Kron-Schritt viele RK-Schritte und werden die Faltungen FFT-basiert organisiert, entsteht ein Vorteil – insbesondere bei großen \(n\). Ferner wird eine komplette Funktion dargestellt.

Bemerkung zur „Weiterdifferenzierbarkeit“

Bei analytischen rechten Seiten führt die Reihenmethode systematisch zu hohen Ableitungsordnungen, weil \[ {}^ky(x_0) = k!\,a_k \] und \(a_k\) durch Faltungen/Kompositionen erzeugt werden. Das ist konzeptionell näher an „Differenzieren durch Rechnen“ (Vektorrechnung/FFT) als an Quadratur mit festem Knotenschema.