, 26 min read

TENDLER: 3. Der Prädiktor

Donald A. Calahan (1972), S.236, schreibt hierzu:

Without prediction $\ldots$ typically 3–5 corrector iterations are required $\ldots$ to satisfy the integration formula $\ldots$ With prediction, convergence is likely in one iteration, although an additional iteration may be necessary to check for convergence. This halving of computational effort should be adequate for incorporating a predictor in an integration algorithm.

Analog äußern sich Byrne/Hindmarsh (1987):

Regardless of the choice of iteration, an initial guess $\bf y_{n(0)}$ is always needed. This is easily obtained by appealing to any of the explicit linear multistep methods.

Der Prädiktor besorgt den für die Korrektoriteration nötigen ersten Startwert. Es wird begründet, warum es lohnenswert ist, speziell gewählte Formeln für den Prädiktor zu verwenden und nicht die Iteration mit anderen, möglicherweise leichter zu berechnenden Werten zu starten.

Man kann leicht beweisen, daß die Prädiktorformeln nicht D-stabil zu sein brauchen, wenn mindestens einmal mit einer D-stabilen Korrektorformel iteriert wird, bzw. das Restverfahren D-stabil ist. Tatsächlich sind die in dem Programm TENDLER verwendeten Prädiktorformeln ab der Ordnung höher oder gleich drei, nicht mehr $D$-stabil. Diese Formeln dürfen daher nicht alleine verwendet werden. Aufgezeigt wird, an welchen Stellen des Programmes TENDLER die Prädiktorformeln eingehen. Zudem wird beschrieben, wie man neue Formeln in das Programm TENDLER einsetzen kann und wie schließlich die Formeln des Prädiktors in dem Programm algorithmisch verwendet wurden. Es zeigt sich dabei, daß die Prädiktor- und Korrektorformeln eng miteinander verwoben sind. Deutlich wird, daß ein Wechseln der Formeln einfach vonstatten gehen kann. Einfache, aber auch nicht-triviale Beispiele werden zur Illustration vorgerechnet. Es ist, anders als in dem Programm STINT, nicht nötig und erforderlich gewisse Besonderheiten eines sehr kompakten Speicherschemas zu beachten und einzuhalten.

Hinzugefügt werden sollte allerdings, daß während der Entwicklungszeit der Programme DIFJMT und STINT, d.h. etwa im Zeitraum von 1970–1978, Speicherersparnisse im Bereich von ca. 1000 Bytes als durchaus lohnenswert erachtet wurden.

Bibliographisch:

  1. Alan C. Hindmarsh (*1942)
  2. George Dennis Byrne (1933–2024), vgl. Nachruf
  3. Byrne/Hindmarsh (1987): Stiff ODE solvers: A review of current and coming attractions, Journal of Comp. Physics, Vol. 70, Issue 1, May 1987
  4. Donald Calahan: “Computer Aided Network Design”, Second Edition, McGraw-Hill, New York St. Louis Düsseldorf $\ldots$, 1972, xvii+350~S.

1. Einbau neuer zyklischer Formeln in das Programm TENDLER

Der Einbau neuer Formeln in das Programm TENDLER ist problemlos möglich. Hat man neue Formeln gefunden, mit denen man praktische Untersuchungen durchführen will, so bietet es sich an, diese in das Programm TENDLER einzusetzen. Hierzu sind zweierlei Dinge zu beachten. Zum einen sind die Formeln selber auszuwechseln und zum anderen die den Formeln zugehörigen Konstanten, wie Fehlerkonstanten und Anzahl der Stufen pro Zyklus.

1. Konstanten. Es sollte allerdings deutlich darauf hingewiesen werden, daß zahlreiche Konstanten, Strategiefaktoren und dergleichen, z.B. beim Konvergenztest, bei der Fehlerkontrolle, bei der Schrittweiten- und Ordnungssteuerung, zwar nicht direkt abhängig sind von der Ordnung und den benutzten Formeln, jedoch das Zusammenwirken aller Einflußgrößen in einem größerem Programm noch zu wenig verstanden ist, als daß man uneingeschränkt sagen könnte, daß diese Konstanten heuristischen Ursprunges gänzlich unabhängig sind von den verwendeten Verfahren. Man wird wohl erwarten, daß sich mit veränderten Formeln auch diese ändern könnten. Wie sich diese nun aber ganz genau ändern, kann im voraus nicht allgemein gesagt werden. I.d.R. ergeben sich diese Konstanten und Strategiefaktoren durch genaues Beobachten von zahlreichen Testläufen mit der verschiedenstartigen Kombination von Steuergrößen und durch Plausibilitätsüberlegungen. Wenn man so will, ergibt sich diese Sammlung von Konstanten, die keine echten “Naturkonstanten” als solche sind, ähnlich wie ein bisher unbekannter Beweis: nicht alleine durch Deduktion aus einem fertigen Axiomensystem.

Um diesen Punkt zu verdeutlichen, ein Beispiel. Der Konvergenztest ist in dem Programm TENDLER gänzlich unabhängig von den Eigenschaften der zugrunde liegenden zyklischen Formeln, jedoch abhängig von der benutzergewählten Genauigkeitsanforderung $\varepsilon$.

In dem Programm STINT ist der Konvergenztest abhängig von den Fehlerfaktoren der Stufen. Für das Programm TENDLER wurde vermeint, daß die Unabhängigkeit des Konvergenztests von den Fehlerfaktoren ein Vorteil ist. Gilt dies nun auch rückwärts für das Programm STINT? Wohl kaum, denn dort wurde gefunden, daß es sich genau umgekehrt verhält, daß also die Aufnahme der Fehlerfaktoren mit in den Konvergenztest günstig ist, siehe Tendler/Bickart/Picel (1978), Seite 350ff. und im Report Tendler/Bickart/Picel (1976), S.25, S.39.

Die Programme STINT und TENDLER verwenden die völlig gleichen zyklischen Formeln, nämlich die sieben steif-stabilen, zyklischen Formeln von J.M. Tendler, siehe Tendler (1973). Allerdings differieren die Schrittweiten- und Ordnungssteuerungen beider Programme und eben die Konvergenztests nicht unerheblich.

Diese Diskrepanz ist nicht gänzlich überraschend. Ähnliche Beobachtungen wurden auch mit den beiden Programmen DIFSUB und GEAR gemacht. Auch hier verwenden die beiden Programme die völlig gleichen Formeln in gleicher Darstellung (Nordsieck-Darstellung), jedoch ist das Programm GEAR das deutlich effizientere.

Wenn also hier von problemloser Integration von neuen Formeln gesprochen wird, dann sind diese Bemerkungen wohl zu beachten. Verschiedene Programme basierend auf ein und derselben Formel, können sich als sehr verschieden erweisen, man vergl. auch die beiden Programme DIFSUB und DE/STEP. Obwohl hier schon die programmiermässigen Unterschiede als solches größer sein müssen, aufgrund der verschiedenen Darstellungsform der Lösungen und deren Ableitungen oder Differenzen, verwenden jedoch beide Programme vom Wesen her die gleiche Formelfamilie.

Ähnliche Überlegungen gelten für DIFBDF vs. DIFSUB, siehe Tendler (1973), und BDF-ODIOUS vs. LSODE, siehe die Dissertation von Peter Tischer von 1983.

  1. Tischer, Peter E.: “The Cyclic Use of Linear Multistep Formulas for the Solution of Stiff Differential Equations”, Ph.D. Thesis, Department of Computer Science, Monash University, Clayton, Victoria, Australia, August 1983, x+180 S.
  2. Tendler, Joel Marvin und Bickart, Theodore A. und Picel, Zdenek: “STINT: STiff ordinary differential equations INTegrator. Program and User Manual”, Electrical and Computer Engineering Department, Syracuse University, Syracuse, New York 13210, USA, and Department of Electrical Engineering, Delft University of Technology, Delft, The Netherlands, Technical Report TR-76-12, December 1976, ii+85 pages

2. Wechselwirkungen. Nun soll an dieser Stelle eine grundsätzliche Problematik beleuchtet werden. Mit der Veränderung gewisser Segmente im Programm TENDLER und einer dann sich anschließenden Überprüfung der Wirksamkeit der Veränderung, kann folgendes Phänomen zu Tage treten. Eine Veränderung, z.B. ein Einbau neuer Formelpaare, von der man aufgrund von Vorüberlegungen eine signifikante Effizienzsteigerung erwartet, kann dann bei einem Test entgegengesetztes Verhalten zeitigen, d.h. Leistungseinbußen aufweisen. Der Schluß, daß die Veränderung, z.B. die neuen Formeln, per se effizienzmindernd wirken, ist jedoch so ohne weitere Überlegungen nicht statthaft. Aufgrund des nicht gänzlich verstandenen Zusammenwirkens der verschiedenen Segmente im Programm TENDLER muß man in die Erwägungen mit einbeziehen, daß eine Veränderung des einen Segmentes, eine Veränderung eines anderen Segmentes notwendig machen könnte. Beispielsweise könnte der Einbau neuer Formeln es viel geeigneter erscheinen lassen, die Schrittweiten- und Ordnungssteuerung grundsätzlich zu verändern. Fände man Formeln, die sagen wir bis zur Ordnung 10 noch $A_\infty^0$-stabil, kleinste Fehlerkonstanten besäßen und noch viel mehr, so wäre der Trend der Schrittweiten- und Ordnungssteuerung zu niedrigen Ordnungen möglicherweise hinderlich. Dieser eben skizzierte Gedankengang bzgl. der Schrittweiten- und Ordnungssteuerung wurde auch tatsächlich so einmal duchprobiert von Gupta (1985) im Programm DSTIFF.

Es wurden die niedrigen Ordnungen nicht mehr so stark bevorzugt. Der erhoffte Erfolg blieb jedoch aus, siehe Gupta (1985).

Die Erläuterungen anhand des Programmes TENDLER haben lediglich illustrativen Charakter. In den Programmen DIFSUB, GEAR, EPISODE, DE/STEP und zahlreichen anderen Programmen ist das Zusammenspiel der verschiedenen Segmente nicht mehr oder minder verstanden.

Es tritt nun jedoch noch eine gravierende Zusatzschwierigkeit hinzu. Glaubt man nun, daß eine Veränderung des einen Segmentes im Programm TENDLER auch eine Veränderung eines anderen Segmentes bedingt, und verfährt man auch entsprechend, d.h. modifiziert man auch weitere Segmente und stellt nun eine Effizienzverbesserung fest, so wäre der Umkehrschluß, daß diese Veränderung nun als solches effizienzanhebend sei, ebenfalls nicht zulässig. Würde man beispielsweise bei einem Einbau neuer Formeln gleichzeitig auch die Fehlerkontrolle und die Schrittweiten- und Ordnungssteuerung ändern, so wäre der Gedanke, daß die neuen Formeln nun günstiger wären, als die zyklischen Formeln von Joel Marvin Tendler, möglicherweise naheliegend, aber dennoch irrig. Hinter alle dem könnte sich letztlich die banale Erkenntnis verbergen, daß ein Programm umso besser wird, desto mehr Pflege man ihm angedeihen läßt.

Denn möglicherweise hätte man das alte Programm in seiner Wirksamkeit ebenso verbessern können und zwar ohne an den Einbau neuer Formeln überhaupt zu denken.

Diese Ausführungen sollen nun nicht den Eindruck erwecken, als ob ein Einbau neuer Formeln als solches schädlich sei.

Die Ausführungen beleuchten lediglich, daß gewisse Folgerungen aus Testläufen, mögen sie auch noch so umfangreich sein, vorsichtig gezogen werden müssen, damit sie nicht in die falsche Richtung weisen. Hinter alle dem steckt das grundsätzliche Dilemma, daß gewisse statistische Tests nicht anwendbar sind, falls Variablen nicht voneinander unabhängig sind.

2. Alternativen zum Prädiktor

1. Um einen groben Überblick zu bekommen, ob es sich überhaupt rechnet einen Prädiktorschritt der Korrektoriteration voranzustellen, seien einige Erläuterungen gegeben. Bekanntermaßen kann selbst ein Korrektor mit leerem Stabilitätsgebiet, mit Hilfe eines vorgeschalteten Prädiktors im Falle eines Picard-$P(EC)^i{E}$-Verfahrens, zu einem Gesamtverfahren werden, welches einen akzeptablen Stabilitätsbereich aufweist. Insbesondere kann durch geschickte Wahl eines Prädiktors, der Stabilitätsbereich eines Picard-$P(EC)^i{E}$-Verfahrens vergrößert werden, siehe Albrecht (1979), §9.7, S.128–131 und Albrecht (1976), §9.5.2.

Diese Eigenschaft ist insoweit besonders hervorhebenswert, da für die $D$-Stabilität eines Picard-$P(EC)^i{E}$-Verfahrens die Prädiktorformel völlig unerheblich ist, insbesondere braucht die Prädiktorformel nicht selber $D$-stabil zu sein. Erkennt man nun das Potential der möglichen Effizienzsteigerung und der Stabilitätsverbesserung durch Benutzung eines Prädiktors, so kann man nach dem besten unter allen Prädiktoren fragen für einen gegebenen Korrektor. Oder aber man fragt nach der besten Kombination von Prädiktor/Korrektor. “Besser” wäre: größeres Stabilitätsgebiet, besseres Rundungsfehlerverhalten, etc. Diese Frage wird hier nirgends beantwortet.

Da sich die Picard Iteration (mit $J=0$) als Spezialfall des modifizierten Newton-Kantorovich Iterationsverfahren auffassen läßt und zusätzlich die Anzahl der Iterationsschritte beim augenblicklich verwendeten Newton-Verfahren im Programm TENDLER sehr gering ist (meistens 1, selten 2, ganz selten 3 und mehr Iterationsschritte), liegt die Vermutung nahe, daß die Ergebnisse, die für Picard-$P(EC)^i{E}$-Verfahren gelten, in wenn auch abgeschwächter Form, ebenso für das modifizierte Newton-Kantorovich-$P(EC)^i{E}$-Verfahren gelten. Erneut taucht die Frage auf, welche Kombination von Prädiktor, Korrektor und Iterationsverfahren unter allen anderen Verfahren das beste und günstigste ist.

“Besser” hier im Sinne von: Minimierung der Anzahl der Iterationen, größere Stabilitsgebiete, Rundungsfehlereigenschaften, u.s.w. Auch diese Frage wird hier nirgends beantwortet.

Nicht selten haben explizite lineare Mehrschrittverfahren erheblich kleinere Stabilitätsgebiete als (von der Ordnung) vergleichbare implizite lineare Mehrschrittverfahren. Dies ist einer der maßgeblichen Gründe für die Verwendung impliziter linearer Verfahren, anstelle rein expliziter Methoden.

Hinzu kommt ferner, daß die Fehlerkonstanten das implizite Verfahren favorisieren. Typisches Beispiel sind hier die Adams-Formeln (Adams/Moulton vers. Adams/Bashforth).

2. Es gibt verschiedene Interpretationsmöglichkeiten für den Sinn der Benutzung eines Prädiktors. Einmal könnte man denken, daß der Prädiktor die Lösung der Differentialgleichung “voraussagen” soll. In diesem Falle könnte man meinen, daß eine explizite Prädiktorformel, eben aufgrund seines kleinen Stabilitätsbereiches (wenn überhaupt) sehr ungenaue Voraussagungen macht. Dem dieser ersten Interpretation mögliche innewohnende Trugschluß, nämlich instabile Stufe gleich ungenaue Näherungenswerte des gesamten mehrstufigen Verfahrens, soll hier nicht weiter betrachtet werden. Ein Picard- oder Newton-$P(EC)^i{E}$-Verfahren ist nämlich als Gesamtverfahren auffaßbar als zusammengesetztes, mehrstufiges Verfahren. Eine andere Deutungsmöglichkeit liegt darin, daß der Prädiktor möglichst nahe der Lösung von $y=h\gamma f(y)+\psi$ liegen sollte.

Wegen der ersten Deutungsmöglichkeit könnte man vorschlagen, statt einer i.a. ohnehin instabilen Prädiktorformel, einfach faute de mieux den zuletzt akzeptierten Näherungswert für die Lösung der Differentialgleichung als Startwert für die Iteration zu wählen. Dieser Vorschlag wird zumeist im Zusammenhang mit dem Newton-$P(EC)^i{E}$-Verfahren gemacht, nicht jedoch im Zusammenhang mit der Picard Iteration. Dennoch, aufgrund der engen Verwandtschaft beider Iterationsvorschriften, scheint es angebracht, diesen Vorschlag des Startwertes anhand der Picard Iteration zu prüfen, wenn man sich der Beschränkung bewußt ist, die damit einhergeht. Hinzu kommt, daß Steifheit und Instabilität von Formeln letztlich aus gleicher Quelle gespeist werden.

Wählt man das explizite Euler-Verfahren $y_{n+1}=y_n+hf_n$ als Prädiktor und das implizite Euler-Verfahren $y_{n+1}=y_n+hf_{n+1}$ als Korrektor, so lieferte für die Differentialgleichung $\dot y=\lambda y$, das $PECE$-Verfahren die Näherungen

$$ y_{n+1}=(1+H+H^2)y_n, $$

mit $H=h\lambda$.

Hingegen entartet der Korrektor mit dem Startwert $y_{n+1}^0=y_n$ zum expliziten Euler-Verfahren, also genau dem Verfahren, dem man durch die veränderte Wahl eines Startwertes zu entrinnen versuchte. Dank des $H^2$-Terms hat das $PECE$-Verfahren die Möglichkeit einen negativen $H$-Wert noch “in's Positive zu ziehen”. Erwartungsgemäß hat das $PECE$-Verfahren einen größeren Stabilitätsbereich und dies obwohl der Prädiktor als Einzelverfahren betrachtet, seinen Stabilitätsbereich schon verlassen haben könnte. Dieses Verhalten ist typisch für $P(EC)^i\{E\}$-Verfahren. Bei Verwendung von $y_{n+1}^0=y_n$ und Einsetzung, hat man ohnehin ein neues explizites Verfahren vor sich, welches mit seiner Konsistenzordnung weit hinter der Konsistenzordnung des Korrektors abfallen könnte und somit zahlreiche Picard Iterationen benötigen könnte.

Auch noch von einer anderen Seite wurde dem Prädiktor Aufmerksamkeit gewidmet. Tischer/Gupta (1983) (s.u.), §2.5.1, S.10–17, und §3.2, S.38–42, untersuchten, ob die Verwendung linear impliziter Prädiktoren eine signifikante Leistungssteigerung in dem Programm ODIOUS bewirken konnte. Auch mit diesem Prädiktor konnte das Programm ODIOUS basierend auf den zyklischen Formeln von Tischer (1983) (Dissertation s.o.) und Tischer/Sacks-Davis (1983) das Programm LSODE in seiner allgemeinen Leistungsfähigkeit nicht klar überbieten. Dies zeigt, w.o. schon erwähnt, nicht die Untauglichkeit dieser Maßnahme als solches. Es deutet jedoch an, daß noch andere Segmente eines Programmes erheblichen Einfluß auf die Effizienz haben und somit eine einseitige Konzentration der Anstrengungen auf den Prädiktor, nicht zum erhofften Ziele führen.

  1. Tischer, Peter E.: “The Cyclic Use of Linear Multistep Formulas for the Solution of Stiff Differential Equations”, Ph.D. Thesis, Department of Computer Science, Monash University, Clayton, Victoria, Australia, August 1983, x+180 S.

  2. Tischer, Peter E. und Gupta, Gopal K.: “A Cyclic Method Stiff ODE Solver”, Technical Report 38, Department of Computer Science, Monash University, Clayton, Victoria, Australia, June 1983, i+45 S.

3. Verwebung von Prädiktor und Korrektor

It may look a bit formidable at first glance; pointers fly most alarmingly. $\ldots$ the doubly linked list is the representational technique of choice in many circumstances.

W.A. Wulf, M. Shaw, L. Flon, P.N. Hilfinger (1981)

1. Anzahl Stufen pro Zyklus. Die Anzahl der Stufen pro Zyklus ist in einem Konstantenvektor abgespeichert, der z.Z. wie folgt aussieht:

stgs_per_order[8] = {*,3,3,3,3,4,4,4};

Die Fehlerfaktoren sind in den beiden Konstantenvektoren abgelegt:

alfa_div_gamma[8][4] = { ... };
errconst[8] = { ... };

Für die Fixpunktiteration wird $\gamma_i=\beta_{ii}/\alpha_{ii}$ benötigt, abgespeichert in beta_div_alfa[8][5].

Für den Hindmarsh-Test und für die modifizierte Newton-Kantorovich Iteration wird meanbetalf[8] benötigt, als Mittelwert der $\gamma_i$. Die hier vorzunehmenden Änderungen hängen allein einzig von den verwendeten Formeln ab und sind offensichtlich.

Benötigen die Formeln zyklusfremde Ableitungen, wie dies beispielsweise bei den Formeln von Tischer (1983) (Dissertation) und Tischer/Sacks-Davis (1983) der Fall ist, bzw. sinngemäß bei den Adams-Formeln, so müssen adapthz() und newbdf() angepasst werden, damit diese Ableitungen für das neue Gitter bereitgestellt werden können.

Weiterhin muß in diesem Falle die Funktion actbdft() modifiziert werden, damit pro erfolgreichem Schritt eine entsprechende Differenzentabelle für die Ableitungen mit aufgebaut wird. Diese Änderungen, die über das reine Ändern von Konstanten hinaus gehen, sind durch entsprechende Kommentare im Programm gekennzeichnet. U.U. ist auch mehr Speicherplatz bereitzustellen.

2. Schreibweise. Nicht mehr ganz so offensichtlich sind die Änderungen für die neuen Formelpaare selber. Hierzu folgen nun Beispiele, wie anhand konkreter Formeln entsprechende Werte gesetzt werden müssen. Es ist zu beachten, daß die einzusetzenden Konstanten vom verwendeten Prädiktor abhängen. Hierzu beobachtet man die folgenden Relationen.

Zur notationellen Vereinfachung sei die Funktion $f$ in autonomer Schreibweise angenommen, und es sei abkürzend geschrieben $z = hf(y).$ Die zugrunde liegende Mehrschrittformel laute

$$ \sum_{i=0}^\kappa \alpha_i y_{n+i-\kappa+1} = h \sum_{i=0}^\kappa \beta_i f_{n+i-\kappa+1} . $$

Mit

$$ y_{n+1} = h\gamma f(y_{n+1})+\psi = \gamma z_{n+1}+\psi $$

erhält man für den Wert $z_{n+1}$ dann

$$ z_{n+1} = {1\over\gamma}(y_{n+1} - \psi), $$

insbesondere

$$ z_{n+1}^0 = {1\over\gamma}(y_{n+1}^0 - \psi). $$

Dies ist die ganz zentrale und entscheidende Formel.

Bei den weiteren Überlegungen wird diese Formel immer wieder auftauchen. Es gilt $z_{n+1}^0$ zu bestimmen, dabei sollen für $y$ rückwärtsgenommene Differenzen verwendet werden.

Der Wert $\psi$ sammelt bei einem linearen Mehrschrittverfahren alle schon vorher berechneten Terme auf, es ist also

$$ \psi = {1\over\alpha_i} \sum_{i=0}^{\kappa-1}\left(-\alpha_iy_{n+i-\kappa+1} + h\beta_ihf_{n+i-\kappa+1}\right). $$

Berücksichtigt man die Stelle der Stufe innerhalb des Zykluses, so hätte man zu schreiben

$$ \psi_{ij} = {1\over\alpha_{ii}} \sum_{j=-\kappa+i}^{i-1} \left( -\alpha_{ij} y_{m\ell+j} + h \beta_{ij} f_{m\ell+j} \right) . $$

Zur Schreibvereinfachung seien Indizes unterdrückt, solange dies nicht zu Mißverständnissen führen kann.

Weiter sei zur Vereinfachung der Schreibweise folgende Vereinbarung eingeführt. Der Term $\psi$, der also die vergangenen Werte der $y_i$ und $z_i$ enthält, wird weiter aufgespalten in die beiden Summanden $\psi_y$ und $\psi_z$, wobei $\psi_y$ die Summe der vergangenen $y_i$ Werte enthält und $\psi_z$ diejenigen Werte für $z_i$.

Z.B. ist für das implizite Euler-Verfahren $\psi_y=y_n$ und $\psi_z=0$.

Es folgen nun eine Reihe von Beispielen, die einem das Gefühl vermitteln, welche Rechnungen durchzuführen sind und wie diese konkret aussehen. Dabei werden auch nebenbei einige praktische Hinweise gegeben.

3. Beispiel für Ordnung 1: Für das implizite Euler-Verfahren

$$ y_{n+1} = y_n + hf_{n+1} = y_n + z_{n+1} $$

mit dem expliziten Euler-Verfahren $y^0_{n+1}=y_n+z_n$ als Prädiktor ergibt sich

$$ z^0_{n+1}=(y^0_{n+1}-\psi)/\gamma=z_n+0y_n, \hbox{ wegen } \psi=y_n . $$

4. Beispiel für Ordnung 2: Bei dem Zyklus zweiter Ordnung, bei dem alle Stufen untereinander gleich sind und zwar gleich der BDF2, gehen die Rechnungen wie folgt vor sich. Die BDF2 lautet

$$ 3y_{n+1}\underbrace{-4y_n+y_{n-1}}_{\displaystyle{-3\psi_y}} = 2z_{n+1} $$

und der Prädiktor sei

$$ y^0_{n+1} = y_{n-1}+2z_n. $$

Es ist also $\gamma=2/3$, und damit rechnet man

$$ \eqalign{ z_{n+1}^0 &= {1\over\gamma}(y^0_{n+1}-\psi)\cr &= {3\over2}\biggl[(y_{n-1}+2z_n)-{1\over3}(4y_n-y_{n-1})\biggr]\cr &= -2y_n+2y_{n-1}+3z_n.\cr } $$

Die führende Stelle ist hier natürlich wieder $p/\gamma$. Der Anteil mit $y$ in $z^0_{n+1}$ wird umgeformt in die Darstellung mit Hilfe rückwärtsgenommener Differenzen durch

$$ \tilde\alpha = \hat P^\top\alpha = \pmatrix{1 & 1\cr 0 & -1\cr}\pmatrix{-2\cr 2\cr} = \pmatrix{0\cr -2\cr}. $$

5. Beispiel für Ordnung 3: Für die erste Stufe im Zyklus der Ordnung $p=3$, also BDF3, mit

$$ 11y_{n+1}\underbrace{-18y_n+9y_{n-1}-2y_{n-2}}_{\displaystyle-11\psi_y} = 6z_{n+1} $$

und dem Prädiktor

$$ y^0_{n+1} = {1\over2}\bigl(-3y_n+6y_{n-1}-y_{n-2}\bigr)+3z_n $$

ergibt sich mit $\psi=\psi_y$ und $p/\gamma=11/2$

$$ \eqalign{ z_{n+1}^0 &= {1\over\gamma}(y_{n+1}^0-\psi)\cr &= {11\over6}\biggl[(-{3\over2}y_n+3y_{n-1}-{1\over2}y_{n-2}) -{1\over11}(18y_n-9y_{n-1}+2y_{n-2})+3z_n\biggr]\cr &= \bigl(-{23\over4}y_n+7y_{n-1}-{5\over4}y_{n-2}\bigr)+{11\over2}z_n.\cr } $$

Schließlich ergibt die Umrechnung in die Darstellung mit Hilfe rückwärtsgenommener Differenzen

$$ \tilde\alpha = \pmatrix{1 & 1 & 1\cr 0 & -1 & -2\cr 0 & 0 & 1\cr} \pmatrix{-23/4\cr 7\cr -5/4\cr} = \pmatrix{0\cr -9/2\cr -5/4\cr}. $$

6. Beispiel für Ordnung 4: Die erste Stufe im Zyklus der Ordnung $p=4$ lautet

$$ 25y_{n+1}\underbrace{-48y_n+36y_{n-1}-16y_{n-2}+3y_{n-3}}_{=-25\psi_y} =12z_{n+1}, $$

welches natürlich gerade die BDF4 ist. Hierbei ist der Wert $\gamma$ definiert durch

$$ \gamma = {\beta_{n+1}\over\alpha_{n+1}}. $$

$\psi_z$ ist Null. In anderer Formulierung läßt sich die obige Formel also schreiben zu

$$ % 25y_{n+1}=12z_{n+1}+25\psi_y\qquad\hbox{also}\qquad y_{n+1} = {12\over25}z_{n+1} + \psi_y. $$

Mit $z_{n+1}^0 = (y_{n+1}^0-\psi)/\gamma$ und dem Prädiktor für $y_{n+1}^0$ zu

$$ y_{n+1}^0 = {1\over3}(-10y_n+18y_{n-1}-6y_{n-2}+y_{n-3})+4z_n, $$

ergibt sich jetzt zunächst für den Anteil mit $z$-Werten der entsprechende Prädiktorwert für $z_{n+1}^0$ zu

$$ z_{n+1}^0\bigl|_z = {1\over\gamma}(y_{n+1}^0\bigl|_z-\psi_z) = {25\over12}(4z_n-0) = {25\over3}z_n, $$

also, wie zu erwarten, der Wert $p/\gamma$, aufgrund der speziellen Wahl des Prädiktors für $y$. Dieser Wert, also $25/3$, wird nun abgespeichert in dem Vektor zpredbeta[3][0][0]. Für den Anteil mit $y$ erhält man mit der gleichen Überlegung

$$ \eqalign{ z_{n+1}^0\bigl|_y &= {1\over\gamma}(y_{n+1}^0\bigl|_y-\psi_y)\cr &= {25\over12}\biggl[{1\over3}(-10y_n+18y_{n-1}-6y_{n-2}+y_{n-3}) -{1\over25}(+48y_n-36y_{n-1}+16y_{n-2}-3y_{n-3})\biggr]\cr &= -{197\over18}y_n+{31\over2}y_{n-1}-{11\over2}y_{n-2} +{17\over18}y_{n-3}.\cr } $$

Dieser Wert muß jetzt noch in die Darstellung der rückwärtsgenommenen Differenzen umgeformt werden.

Zwischen der Darstellung in Ordinatendarstellung und Darstellung in Form rückwärtsgenommener Differenzen herrscht der Zusammenhang $Y=\hat P\nabla Y$, oder ausgeschrieben

$$ \pmatrix{y_n\cr y_{n-1}\cr y_{n-2}\cr \vdots\cr y_0\cr} = \pmatrix{ 1 & & & & \cr 1 & -1 & & & 0\cr 1 & -2 & 1 & & \cr \vdots & \vdots & \vdots & \ddots & \cr 1 & -n & n\choose2 & \ldots & (-1)^n\cr} \pmatrix{\nabla^0y_n\cr \nabla^1y_n\cr \nabla^2y_n\cr \vdots\cr \nabla^ny_n\cr} . $$

Mit dieser Beziehung berechnet man nun, mit den Vektoren

$$ Y_n := \pmatrix{y_n\cr y_{n-1}\cr y_{n-2}\cr y_{n-3}\cr} \qquad\hbox{und}\qquad \alpha := \pmatrix{-197/18\cr 31/2\cr -11/2\cr 17/18\cr} $$

und der Gleichungskette

$$ z_{n+1}^0 = \alpha^\top Y_n = (\hat P^\top\alpha)^\top \nabla Y_n = \tilde\alpha^\top \nabla Y_n, $$

schließlich den Wert $\tilde\alpha$ einfach durch $\tilde\alpha=\hat P^\top\alpha$, also

$$ \pmatrix{ 1 & 1 & 1 & 1\cr 0 & -1 & -2 & -3\cr 0 & 0 & 1 & 3\cr 0 & 0 & 0 & -1\cr} \pmatrix{-197/18\cr 31/2\cr -11/2\cr 17/18\cr} = \pmatrix{0\cr -22/3\cr -8/3\cr -17/18\cr}. $$

Dieser Wert wird dann auch tatsächlich abgespeichert. Die einzige Erklärung, die jetzt noch verbleibt zu geben, ist, daß aufgrund der Konsistenz des Prädiktors mindestens der Ordnung $p=1$, die oberste Komponente immer Null sein muß, wegen der speziellen Beschaffenheit der Matrix $\hat P^\top$, nämlich der zuobersten “Einerzeile” wie bei der Konsistenzmatrix $C_{p,k}$. Mit den offensichtlichen Vereinfachungen, werden diese Rechnungen am besten mit Formelmanipulationssystemen durchgeführt.

In der Sprache muMATH z.B. würde man eingeben

AS: P.(Y-PSIY)/GAMMA;

7. Beispiel für Ordnung 6: Die zweite Stufe im Zyklus der Ordnung $p=6$ lautet

$$ 2930y_{n+1}\underbrace{-7277y_n+9150y_{n-1}-8100y_{n-2}+4550y_{n-3} -1455y_{n-4}+202y_{n-5}}_{\displaystyle -2930\psi_y} = 1200z_{n+1}-60z_n. $$

Der Prädiktor für diese Stufe mit entsprechender Konsistenzordnung ist

$$ y^0_{n+1} = {1\over10}(-77y_n+150y_{n-1}-100y_{n-2}+50y_{n-3}-15y_{n-4} +2y_{n-5})+6z_n $$

Es ist also $\gamma=1200/2930$ und $\psi_z=-60/2930\cdot z_n$. Der Anteil mit $z$ für den Prädiktor $z^0_{n+1}$ ist somit

$$ z_{n+1}^0\bigl|_z = {1\over\gamma}(y^0_{n+1}\bigl|_z-\psi_z) = {2930\over1200}(6z_n+{60\over2930}z_n) = {147\over10}z_n $$

und für den Anteil mit $y$ erhält man

$$ z_{n+1}^0\bigl|_y = {1\over\gamma}(y_{n+1}^0\bigl|_y-\psi_y) = -{4973\over200}y_n+{177\over4}y_{n-1}-{187\over6}y_{n-2} +16y_{n-3}-{117\over24}y_{n-4}+{197\over300}y_{n-5}. $$

Umrechnung dieses Wertes in Form von rückwärtsgenommenen Differenzen liefert

$$ \tilde\alpha = \pmatrix{ 1 & 1 & 1 & 1 & 1 & 1\cr 0 & -1 & -2 & -3 & -4 & -5\cr 0 & 0 & 1 & 3 & 6 & 10\cr 0 & 0 & 0 & -1 & -4 & -10\cr 0 & 0 & 0 & 0 & 1 & 5\cr 0 & 0 & 0 & 0 & 0 & -1\cr} \pmatrix{-4973/200\cr 177/4\cr -187/6\cr 16\cr -117/24\cr 197/300\cr} = \pmatrix{0\cr -137/10\cr -117/20\cr -46/15\cr -191/120\cr -197/300\cr}. $$

8. Zyklusfremde Näherungen. Bei den Formeln von Tendler, siehe auch Tendler (1973), sind die Formeln für $p\le 3$ stets von der Gestalt

$$ A_0u_n+A_1u_{n+1} = B_0z_n+B_1z_{n+1}, \qquad\hbox{mit}\quad B_0=0, $$

und für $p\ge 4$:

$$ A_0u_{n-1}+A_1u_n+A_2u_{n+1} = B_0z_{n-1}+B_1z_n+B_2z_{n+1}, \qquad\hbox{mit}\quad B_0=B_1=0. $$

Also die zyklischen Formeln benutzen keine zyklusfremden Näherungen für Ableitungen. Die verwendeten Ableitungen stammen sämtlich aus dem aktuellen Zyklus. Die Prädiktoren für alle Ordnungen benutzen nur eine einzige Ableitung aus dem vorhergehenden Zyklus und zwar die letztmögliche Ableitungsnäherung. Diese spezielle Wahl der Matrizen $B_i$ hat nun entscheidende Konsequenzen: Ist das Verfahren $A[\alpha]$- bzw. $S[\delta]$-stabil, so ist es automatisch auch $A_\infty^0[\alpha]$- bzw. $S_\infty^0[\delta]$-stabil.

I.d.S. sind also die Formeln von J.M. Tendler die natürliche Verallgemeinerung der BDF. Dieser Automatismus ist ein äußerst einfacher Weg zur Sicherstellung von $A_\infty^0[\alpha]$- bzw. $S_\infty^0[\delta]$-Stabilität.

Von mindestens ebenso großen Nutzen ist der enorme rechnerische, speicherplatzmässige und programmiermässige Vorteil. Die innerhalb eines Zykluses generierten Näherungswerte für die Ableitungen können nach Beendigung des Zykluses einfach vergessen werden, also überschrieben werden und vielmehr: Sie bedürfen keinerlei Aktualisierungsaufwand, wie beispielsweise die rückwärtsgenommenen Differenzen der $y$-Werte. Jede neuberechnete Stufennäherung innerhalb eines Zykluses fügt einen weiteren $y$-Wert hinzu und daher muß die Differenzentabelle im Speicher um eine Spalte erweitert werden. Dies sind die Werte $\nabla y, \nabla^2y, \ldots, \nabla^py$. Dieser Rechenaufwand entfällt völlig für die Ableitungsnäherungen bei den Formeln von Tendler (1973), insbesondere bei den BDF, wenn man sie in der Darstellung rückwärtsgenommener Differenzen o.ä. programmiert.

Von weiterem großen Vorteil ist aber nicht nur der verringerte Aktualisierungsaufwand pro Schritt, sondern der halbierte Aufwand bei der Interpolation von Zwischengitterpunkten bei einem Schrittweitenwechsel. Benötigt nämlich die Formel zyklusfremde Ableitungswerte, so sind natürlich auch diese auf einem neuen Gitter zur Verfügung zu stellen. Genau dies ist eine der Aufgaben der Differenzentabelle. Man könnte zwar auch daran denken ohne Differenzentabelle zu interpolieren, jedoch ist eine mehrfache Interpolation erforderlich, sodaß sich eine Differenzentabelle besonders anbietet. Diese Differenzentabelle für die Ableitungen belegt aber nicht nur zusätzlichen Speicherplatz, benötigt zusätzliche Rechenzeit, sondern läßt sich auch nicht so einfach vermeiden, da bei einer Schrittweitenverkleinerung auf jeden Fall interpoliert werden muß.

Eine Schrittweitenverkleinerung kann nach jeder Stufe vonnöten sein, je nach dem wie man die Fehlerkontrolle organisiert hat. Die Summe dieser Nachteile, falls also $B_0\ne0$ bzw. $\left|B_0\right|^2+\left|B_1\right|^2\ne0$, kann nur dann aufgewogen werden, wenn die Stabilitäts- und Fehlereigenschaften neuer Formeln signifikant besser sind, als diejenigen von Tendler (1973). Die Formeln von Tischer (1983) (Dissertation) und Tischer/Sacks-Davis (1983) sind auch tatsächlich keine Verallgemeinerung der BDF. Bei diesen zweistufigen zyklischen Verfahren sind die $B_i\ne0$ ($i$ geeignet).

4. Speicherung des Prädiktors im Programm TENDLER

Pointers have been lumped with the goto statement as a marvelous way to create impossible-to-understand programs. This is certainly true when they are used carelessly, and it is easy to create pointers that point somewhere unexpected. With discipline, however, pointers can also be used to achieve clarity and simplicity. This is the aspect that we will try to illustrate.

B.W. Kernighan, D.M. Ritchie (1978)

Da für den nicht-steifen Falle, also im Programm TENDLER bei Picard Iteration, eine Prädiktorformel sinnvoll ist, wurde dieses Vorgehen sinngemäß übertragen für den steifen Fall, also im Programm TENDLER Verwendung eines modifizierten Newton-Kantorovich Iterationsverfahren. Also, sowohl im steifen als auch im nicht-steifen Modus werden die gleichen Prädiktorformeln verwendet.

1. Prädiktor für $y$. Die Variable ypred[][] ist als static real ypred[7][7] definiert. Vorbesetzt wird sie mit den Werten der expliziten BDF. Diese Werte werden niemals im Programm verändert. Im Sine von C++ sind sie also static const real.

$$ \begin{array}{cccccc} \displaystyle 1 & \cr \displaystyle 1 & \displaystyle{(1-2)\over1}\cr \displaystyle 1 & \displaystyle{(1-3)\over1} & \displaystyle{(2-3)\over2}\cr \displaystyle 1 & \displaystyle{(1-4)\over1} & \displaystyle{(2-4)\over2} & \displaystyle{(3-4)\over3}\cr \displaystyle 1 & \displaystyle{(1-5)\over1} & \displaystyle{(2-5)\over2} & \displaystyle{(3-5)\over3} & \displaystyle{(4-5)\over4}\cr \displaystyle 1 & \displaystyle{(1-6)\over1} & \displaystyle{(2-6)\over2} & \displaystyle{(3-6)\over3} & \displaystyle{(4-6)\over4} & \displaystyle{(5-6)\over5}\cr \displaystyle 1 & \displaystyle{(1-7)\over1} & \displaystyle{(2-7)\over2} & \displaystyle{(3-7)\over3} & \displaystyle{(4-7)\over4} & \displaystyle{(5-7)\over5} & \displaystyle{(6-7)\over6}\cr \end{array} $$

Für die Initialisierung der Konstanten zpredalfa[][][] und zpredbeta[][][] ist folgendes zu beachten. Die Einträge liegen dabei “linksbündig” vor, da sie ja später bei einer Matrix-Vektor-Multiplikation mit der Matrix

$$ \left(y,{\mskip 3mu} \nabla y,{\mskip 3mu} \nabla^2 y,{\mskip 3mu} \ldots, \nabla^{p-1} y\right) $$

gebraucht werden.

Es sei daran erinnert, daß die rückwärtsgenommene Differenz $(p-1)$-ter Ordnung insgesamt $p$ Werte benutzt.

Die Bedeutung der Indizes für $\mathtt{zpredalfa[}p\mathtt{][}\ell\mathtt{][}i\mathtt{]}$ ist nun:

  1. $p$ für die Ordnung des augenblicklich benutzten Verfahrens,
  2. $\ell$ für die Stufe und
  3. $i$ für den $i$-ten Koeffizienten des Verfahrens.

2. Prädiktor für $z$. Die Variable zpredalfa[][][] ist definiert als static real zpredalfa[7][4][6] und verhält sich in ihrem Zugriff wie die Variable ypred[][].

Ein Stern * deutet an, daß hier ein Wert gesetzt wird, der als unmöglich erkannt wird, falls diese Speicherstelle dennoch angesprochen wird, zumindestens im DEBUG Modus.

Der Wert für $\nabla^0y$ braucht natürlich nicht abgespeichert werden, da er aufgrund der Konsistenz des Prädiktors verschwindet.

Ordnung 1: Für zpredalfa[0][][] wird nun die folgende konstante Matrix fest abgelegt:

$$ \begin{array}{cccccc} * & * & * & * & * & * \cr % p = 1 * & * & * & * & * & * \cr * & * & * & * & * & * \cr * & * & * & * & * & * \cr \end{array} $$

Ordnung 2: Die Matrix zpredalfa[1][][] enthält die Teilmatrix:

$$ \begin{array}{cccccc} -2 & * & * & * & * & * \cr % p = 2 -2 & * & * & * & * & * \cr -2 & * & * & * & * & * \cr * & * & * & * & * & * \cr \end{array} $$

Ordnung 3: Für zpredalfa[2][][] werden die Konstanten wie folgt gesetzt:

$$ \begin{array}{cccccc} \displaystyle{-{9\over2}} & \displaystyle{-{5\over4}} & * & * & * & * \cr % p = 3 \displaystyle{-{9\over2}} & \displaystyle{-{5\over4}} & * & * & * & * \cr \displaystyle{-{15\over2}} & \displaystyle{-{3\over4}} & * & * & * & * \cr * & * & * & * & * & * \cr \end{array} $$

Ordnung 4: Für zpredalfa[3][][] die Werte:

$$ \begin{array}{cccccc} \displaystyle{-{22\over3}} & \displaystyle{-{8\over3}} & \displaystyle{-{17\over18}} & * & * & * \cr % p = 4 \displaystyle{-{22\over3}} & \displaystyle{-{8\over3}} & \displaystyle{-{17\over18}} & * & * & * \cr \displaystyle{-9} & \displaystyle{-{9\over4}} & \displaystyle{-{7\over8}} & * & * & * \cr * & * & * & * & * & * \cr \end{array} $$

Ordnung 5: zpredalfa[4][][]

$$ \begin{array}{cccccc} \displaystyle{-{125\over12}} & \displaystyle{-{101\over24}} & \displaystyle{-{71\over36}} & \displaystyle{-{37\over48}} & * & * \cr % p = 5 \displaystyle{-{125\over12}} & \displaystyle{-{101\over24}} & \displaystyle{-{71\over36}} & \displaystyle{-{37\over48}} & * & * \cr \displaystyle{-{253\over24}} & \displaystyle{-{1001\over240}} & \displaystyle{-{707\over360}} & \displaystyle{-{123\over160}} & * & * \cr \displaystyle{-{57\over4}} & \displaystyle{-{25\over8}} & \displaystyle{-{17\over10}} & \displaystyle{-{169\over240}} & * & * \cr \end{array} $$

Ordnung 6: Für zpredalfa[5][][] werden besetzt:

$$ \begin{array}{cccccc} \displaystyle{-{137\over10}} & \displaystyle{-{117\over20}} & \displaystyle{-{46\over15}} & \displaystyle{-{191\over120}} & \displaystyle{-{197\over300}} & * \cr % p = 6 \displaystyle{-{137\over10}} & \displaystyle{-{117\over20}} & \displaystyle{-{46\over15}} & \displaystyle{-{191\over120}} & \displaystyle{-{197\over300}} & * \cr % \lt ------- \displaystyle{-{353\over25}} & \displaystyle{-{571\over100}} & \displaystyle{-{1819\over600}} & \displaystyle{-{79\over50}} & \displaystyle{-{3919\over6000}} & * \cr \displaystyle{-{3529\over200}} & \displaystyle{-{1889\over400}} & \displaystyle{-{1609\over600}} & \displaystyle{-{3527\over2400}} & \displaystyle{-{931\over1500}} & * \cr \end{array} $$

Ordnung 7: Schließlich für zpredalfa[6][][] werden die folgenden Werte in die Teilmatrix hineingeschrieben.

$$ \begin{array}{cccccc} \displaystyle{-{343\over20}} & \displaystyle{-{303\over40}} & \displaystyle{-{253\over60}} & \displaystyle{-{589\over240}} & \displaystyle{-{101\over75}} & \displaystyle{-{23\over40}} \cr % p = 7 \displaystyle{-{343\over20}} & \displaystyle{-{303\over40}} & \displaystyle{-{253\over60}} & \displaystyle{-{589\over240}} & \displaystyle{-{101\over75}} & \displaystyle{-{23\over40}} \cr \displaystyle{-{266\over15}} & \displaystyle{-{221\over30}} & \displaystyle{-{749\over180}} & \displaystyle{-{73\over30}} & \displaystyle{-{803\over600}} & \displaystyle{-{103\over180}} \cr \displaystyle{-{1316\over75}} & \displaystyle{-{1151\over150}} & \displaystyle{-{3689\over900}} & \displaystyle{-{121\over50}} & \displaystyle{-{4003\over3000}} & \displaystyle{-{257\over450}} \cr \end{array} $$

3. Es folgt nun die Variable zpredbeta[][][] definiert zu

static real zpredbeta[7][4][4];

Die Reihenfolge der Indizes ist wie bei der Variablen zpredalfa[][][].

Die Variable zpredbeta[][][] ist wie folgt initialisiert.

Ordnung 1: zpredbeta[0][][]

$$ \begin{array}{cccc} 1 & * & * & * \cr % p = 1 0 & 1 & * & * \cr 0 & 0 & 1 & * \cr * & * & * & * \cr \end{array} $$

Ordnung 2: zpredbeta[1][][]

$$ \begin{array}{cccc} 3 & * & * & * \cr % p = 2 0 & 3 & * & * \cr 0 & 0 & 3 & * \cr * & * & * & * \cr \end{array} $$

Ordnung 3: zpredbeta[2][][]

$$ \begin{array}{cccc} \displaystyle{11\over2} & * & * & * \cr % p = 3 0 & \displaystyle{11\over2} & * & * \cr 0 & 2 & \displaystyle{13\over2} & * \cr * & * & * & * \cr \end{array} $$

Ordnung 4: zpredbeta[3][][]

$$ \begin{array}{cccc} \displaystyle{25\over3} & * & * & * \cr % p = 4 0 & \displaystyle{25\over3} & * & * \cr 0 & \displaystyle{5\over4} & \displaystyle{35\over4} & * \cr * & * & * & * \cr \end{array} $$

Ordnung 5: zpredbeta[4][][]

$$ \begin{array}{cccc} \displaystyle{137\over12} & * & * & * \cr % p = 5 0 & \displaystyle{137\over12} & * & * \cr 0 & \displaystyle{1\over10} & \displaystyle{1373\over120} &* \cr 0 & \displaystyle{-{1\over20}} & 3.1 & \displaystyle{61\over5}\cr \end{array} $$

Ordnung 6 zpredbeta[5][][]

$$ \begin{array}{cccc} \displaystyle{147\over10} & * & * & * \cr % p = 6 0 & \displaystyle{147\over10} & * & * \cr 0 & \displaystyle{7\over20} & \displaystyle{1477\over100} & * \cr 0 & \displaystyle{-{3\over20}} & \displaystyle{17\over5} & \displaystyle{3079\over200}\cr \end{array} $$

Ordnung 7: zpredbeta[6][][]

$$ \begin{array}{cccc} \displaystyle{363\over20} & * & * & * \cr % p = 7 0 & \displaystyle{363\over20} & * & * \cr 0 & \displaystyle{1\over2} & \displaystyle{547\over30} & * \cr 0 & \displaystyle{-{1\over5}} & \displaystyle{1\over2} & \displaystyle{2737\over150}\cr \end{array} $$