, 25 min read

Die verwendeten zyklischen Formeln im Programm TENDLER

Original post is here eklausmeier.goip.de/blog/2025/02-24-die-verwendeten-zyklischen-formeln-im-programm-tendler.


It is entirely possible, and even very likely, in view of the fact that the parameters Cij were “manually” set in an interactive computer search for the maximum α, that better and simpler cyclic composite methods exist. An interested reader is invited to try the development of such methods on his own. For those, who are satisfied with the methods obtained thus far, in what follows, we describe how they are implemented in a numerical algorithm.

J.M. Tendler, Th.A. Bickart, Z. Picel (1976)

In diesem Abschnitt werden sämtliche benutzten Formeln vorgestellt, die in dem Programm TENDLER Verwendung gefunden haben. Zuerst kommen die Korrektorformeln und anschliessend die Prädiktorformeln.

Angegeben werden auch mehrere Fehlerfaktoren und Fehlerkonstanten und natürlich auch die Werte, die maßgeblich die Stabilitätseigenschaften quantifizieren, also die Werte α und δ bei der A[α]-Stabilität und S[δ]-Stabilität. Hingewiesen wird auch auf Charakteristika und auf den prinzipiellen Aufbau der zyklischen Formeln. Diskutiert wird ebenfalls, welche Bedingungen Stabilität bei unendlich in einfachster Weise garantieren.

Mit angeführt werden auch mehrere Darstellungen für diese Formelpaare; nämlich einmal in der gewöhnlichen Darstellung und das andere Mal in Form von rückwärtsgenommenen Differenzen. Die letztgenannte Form ist auch diejenige Form, die auch in dem Programm TENDLER verwendet wird. Die Darstellung in Form von rückwärtsgenommenen Differenzen bietet Vorteile insbesondere bei der Fehlerkontrolle und bei der Schrittweiten- und Ordnungssteuerung. Dort stellen sie gerade Näherungen für entsprechend benötigte Ableitungswerte dar. Die Korrektoriteration hingegen verhält sich neutral gegenüber der Darstellungsart. Ebenso der Prädiktor. In beiden Darstellungsarten, also sowohl der Darstellung in Form von rückwärtsgenommenen Differenzen als auch in gewöhnlicher Ordinatenform, sind die Prädiktorrechnungen einfach durchführbar. Die Prädiktorrechnungen haben einen häufig unterschätzten Einfluß auf die Gesamtrechenzeit.

1. Die sieben zyklischen Korrektorformeln von J.M. Tendler

At first it was hoped that high order A-stable cyclic methods could be found. This goal was not achieved. However, the new integration methods are applicable to the solution of stiff systems. Moreover they exhibit better stability properties than the backward differentiation formulas in Gear's widely accepted algorithm.

J.M. Tendler (1973)

1. Die Verfahren der Ordnung p=1,,7, die J.M. Tendler_{Tendler, Joel Marvin} im Rahmen seiner Dissertation 1973 neu entwickelt hat, werden hier nun angegeben.

Für die an anderer Stelle genauer erläuterte Fehlerkonstante von Henrici

C=vγvρ(1)w,{vρ(1)=0,v0,ρ(1)w=0,w0,

werden auch gleichzeitig die Linkseigenvektoren vp, für p=1,,7 der entsprechenden Matrizen aufgeführt. Der Rechtseigenvektor ist natürlich immer wp=(1,,1), für p=1,,7. Zu dieser Konstanten vgl. man auch den Aufsatz von Albrecht (1985).

Die beiden ersten Verfahren, also die der Ordnung 1 und 2, sind nichts anderes als die zyklische Wiederholung ein und der selben Formel, hier also des impliziten Euler-Verfahrens für die Ordnung 1 und BDF2 für das Verfahren der Ordnung 2.

Bibliographisch:

  1. Joel Marvin Tendler: “A Stiffly Stable Integration Process Using Cyclic Composite Methods”, Ph.D. Diss., Syracuse University, Syracuse, New York, 26-Feb-1973
  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, {\it ii/}+85~S.
  3. Theodore A. Bickart (1936—2023), obituary
  4. Peter Albrecht, 1985: “Numerical Treatment of O.D.E.s: The Theory of A-Methods”
  5. Euler, Leonhard (1707—1783)
  6. Henrici, Peter Karl Eugen (1923—1987)
  7. Zdenek Picel

2. Die Formeln der Ordnung p=1,,7 sind stets von der Gestalt

j=k+1αijym+j=hj=k+1βijfm+j,i=1,,.

In Matrixschreibweise erhält man z.B. für =3 und k=2, also für ein dreistufiges Verfahren mit 2 Startwerten, ausgeschrieben:

(α1,1α1,0α1,1α1,2α1,3β1,1β1,0β1,1β1,2β1,3α2,1α2,0α2,1α2,2α2,3β2,1β2,0β2,1β2,2β2,3α3,1α3,0α3,1α3,2α3,3β3,1β3,0β3,1β3,2β3,3)(y3m1y3my3m+1y3m+2y3m+3hf3m1hf3mhf3m+1hf3m+2hf3m+3)=(000).

Aus Übersichtlichkeitsgründen werden die sieben nachfolgenden zyklischen Verfahren in der folgenden transponierten Schreibweise angegeben.

Später wird mit den Tableaux in transponierter Darstellung auch in dieser Form gerechnet.

pi=1,,jαijjβij

3. Für die zyklische, dreifache Wiederholung des impliziten Euler-Verfahrens

ym+1=ym+hfm+1

und für die zyklische, dreifache Wiederholung der BDF2

3ym+14ym+ym1=2hfm+1,

ergeben sich die beiden transponierten Tableaux für die Verfahren wie untenstehend.

Ordnung 1.

p=112301001110201130010000110020103001

Ordnung 2.

p=21231100041013412034300310000000120020203002

Für die Links- und Rechtseigenvektoren v1, v2, w1 und w2 erhält man

v1=(111),v2=(111),w1=w2=(111)

und für die Fehlerkonstanten von Henrici ergeben sich damit

C1=32,C2=1.

Diese beiden Verfahren werden sowohl in dem Programm STINT als auch in dem neu entwickelten Programm TENDLER in dieser Form benutzt. Da beide Verfahren, also sowohl das implizite Euler-Verfahren und auch BDF2, optimale, also nicht verbesserbare Stabilitätseigenschaften A0[α] und S0[δ] besitzen, wurden sie unverändert so übernommen.

4. Für p=k=3, p=k=4 und p=k=5 konnte Tendler die folgenden Werte finden. Bei allen drei Verfahren sind wieder die beiden ersten Stufen der Zyklen BDF3, BDF4, bzw. BDF5. Nur die letzte Stufe sorgt, bzw. die letzten beiden Stufen sorgen für die Vergrößerung des Widlund-Winkels α bei der A[α]-Stabilität und der Verkleinerung der Widlund-Distanz δ bei der S[δ]-Stabilität.

Ordnung 3.

p=312322001920018901111892011123003200010000000160420643002

Ordnung 4.

p=41233300216301361611048364812548216202527230093300020001000000011206020124830048

Ordnung 5.

p=51234412000375120022007511801300200735133030030019407801137300298016802013730305470300137355954000115840000300002000010000000001600603020600186030060015304000600

Die Vektoren vi mit viρ(1)=0, die Vektoren wi mit ρ(1)wi=0 und die Fehlerkonstanten von Henrici sind

v3=(579),v4=(12112521),w3=w4=(111),C3=154,C4=667470.

Wäre beim Verfahren dritter Ordnung auch noch die letzte Stufe gleich BDF3, also das gesamte Verfahren dreimal BDF3 hintereinander, so erhielte man für die Fehlerkonstante natürlich C~3=3/4. Der Linkseigenvektor von ρ(1), der zugehörige Rechtseigenvektor und die Fehlerkonstante von Henrici zu dem Verfahren 5.ter Ordnung berechnen sich zu

v5=(53232937719518548547211142),w5=(1111),C5=v5γv5ρ(1)w5=104982866620040151.6931623

5. Man beachte bei allen diesen obigen Formeln, als auch bei den noch folgenden Formeln, den grundlegenden Aufbau. In gewisser Hinsicht wird das Potential, welches man prinzipiell durch Einführung der zyklischen Verfahren erhält, nicht in voller Gänze ausgeschöpft, jedoch stellt sich die vermeintliche Einschränkung unter anderen Gesichtspunkten (programmiermässige Einfachheit, Rechenbedarf, Speicherplatzbedarf) als gar keine Einschränkung heraus, sondern als echter Vorteil. Stets sind die Formeln von der Gestalt,

ym+i000000zm+i=hfm+i

D.h. es handelt sich immer um die zyklische Hintereinanderschaltung linearer Mehrschrittverfahren mit fester Anzahl von Startwerten. Insbesondere heißt dies z.B. für die letzte Stufe, daß nicht alle möglichen k+(1) Werte benutzt werden, sondern tatsächlich nur k, oder bei der Ordnung p=3 sogar nur k1. Erneut wird deutlich, daß es sich hier um eine echte Verallgemeinerung der BDF auf Mehrstufigkeit handelt. Insbesondere werden also nicht alle möglichen zurückliegenden Ableitungswerte verwendet. Aufgrund des Einbein-Charakters der BDF sind natürlich auch andere Verallgemeinerungen denkbar. _{Tischer, Peter E.}Tischer (1983)2 und {Tischer, Peter E.}{Sacks-Davis, Ron}Tischer/Sacks-Davis (1983)3 verfahren hier etwas anders. Hier werden in der letzten Stufe der Zyklen ein Startwert zusätzlich verwendet. Neben den potentiell besseren Stabilitätseigenschaften von zyklischen Verfahren ist dies ein weiterer, möglicher Vorteil von zyklischen Verfahren gegenüber einstufigen linearen Mehrschrittverfahren.

Beachtenswert ist auch die Benutzung der skalierten Ableitungswerte hfm+i. Wie bei den ym+i werden hier nicht alle möglichen Werte wirklich benutzt. Zum maßgeblichen Teil hängt dies an der A0[α]- bzw. S0[δ]-Stabilität. Dennoch muß man hierbei im Auge behalten, daß man diese zusätzlichen Werte auch zu speichern hat. Insbesondere hat man auch für den Fall vorzusorgen, daß ein gesamter Zyklus zurückgewiesen wird und somit noch weitere Werte abgespeichert werden müssen. Darüber hinaus entstehen Rechnungen für die skalierten Ableitungen bei der Interpolation für evtl. benötigte Zwischengitterpunkte bei einem Schrittweitenwechsel.

6. Für die Ordnung p=k=6, wo nur die erste Stufe des Zykluses gleich der BDF6 ist und schließlich für p=k=7, für die die BDF noch niemals D-stabil ist, konnte Tendler die folgenden zyklischen, steif-stabilen Verfahren finden. Bei p=k=7 sind die ersten beiden Stufen wieder die BDF7, wobei man dann auch erneut ein Beispiel dafür sieht, daß die ^{zyklische Kombination klassischer (instabiler) Mehrschrittverfahren} das Gesamtverhalten gegenüber dem Einzelverhalten u.U. verbessern kann.

Ordnung 6.

p=6123451000047220200322514551950240045501399285145081004340203903609150754062251147727789051036020293074451845530029441486540002299500004000030000200001000000000160604201802012006040803001200468040001200

Ordnung 7.

p=7123466000054906000417644902100336751764172277424900367562356349144104900131002298802940441017650481601108929401771066290201089112976815930028604236440009748600005000040000300002000010000000001420060084020420186021003001200840040004200

Die entsprechenden Eigenvektoren zu ρ(1) und die Fehlerkonstante von Henrici sind hierbei

v6=(162594313955842346091601830530),w6=(1111),C6=v6γv6ρ(1)w6=21342463130769311.6320697.

Für das letzte siebente Verfahren lauten die Eigenvektoren und die Fehlerkonstante von Henrici

v7=(4749873033327441999731007530),w7=(1111),C7=v7γv7ρ(1)w7=85572910112500181750.684573.

2. Bemerkungen zu den Formeln von J.M. Tendler

1. Man könnte es als entscheidenden Nachteil ansehen, daß die drei Verfahren der Ordnungen p=5, p=6 und p=7 nun 4 Stufen verlangen. Umfangreiche Testläufe mit dem Programm TENDLER zeigten nun jedoch, daß bei einer geringfügigen Modifikation der Schrittweiten- und Ordnungssteuerung Schrittzurückweisungen innerhalb, oder am Ende, eines Zykluses “fast nie” auftreten. Zudem sind Zurückweisungen in dem Programm TENDLER ohnehin die Ausnahme. Die Testläufe legten sogar nahe noch weiter zu gehen. Zyklen mit 5 Stufen wären völlig unproblematisch in das Programm TENDLER zu integrieren. Diese Erfahrungen stehen im Gegensatz zu den _Entwurfs_zielen bei den zyklischen Verfahren von _{Tischer, Peter E.}Tischer (1983)2 und _{Sacks-Davis, Ron}Tischer/Sacks-Davis (1983)3, nicht jedoch zu den praktischen Erfahrungen, die _{Tischer, Peter E.}Tischer (1983)2, {Tischer, Peter E.}{Gupta, Gopal K.}Tischer/Gupta (1983)1 und {Tischer, Peter E.}{Gupta, Gopal K.}Tischer/Gupta (1985)2 mit dem Programm ODIOUS gemacht haben. %

2. Auffällig ist bei den αij-Werten die ^{Vorzeichenalternierung entlang herunter einer Spalte} und bei den βij-Werten sticht ins Auge, die “Konzentration” der Minus-Zeichen “in der Mitte”, falls nicht dort Nullen stehen, wiederum entlang herunter einer Spalte.

Stets sind die ersten beiden Stufe nichts anderes als die BDF, außer bei dem Verfahren sechster Ordnung.

Für die Abhängigkeit der Stufenzahl von der Ordnung p erhält man übersichtlich zusammengefasst,

p12345673333444

welche man auch mit der Formel =3+p/5 berechnen könnte (1p7).

3. Für alle Formeln von _{Tendler, Joel Marvin}Tendler (1973) gilt, daß zi=hfi für im nicht benötigt werden, man vergl. hier auch das vorher angegebene Sternchen-Schema. Stets ist also B0=0, bei invertierbarer Matrix B1. Hiermit sichert man automatisch, daß die Nullstellen des charakteristischen Polynoms detQ(μ,H) für ReH gegen Null streben, wegen

detQ(μ,H)=det[(A1μ+A0)H(B1μ+B0)]ReHdet(B1μ+B0),

(in nicht ganz korrekter Schreibweise) und somit aufgrund von B0=0 dann

detQ(μ,)=±μsdetB1.

Damit enthält das Spektrum des Matrixpolynoms Q(μ,) nur die Null, falls B1 invertierbar vorausgesetzt wird.

Dies ist nicht der einzige mögliche Weg sicherzustellen, daß die Eigenwerte bei ReH= verschwinden. Wichtig ist lediglich, daß das Spektrum des Matrixpolynoms

i=0κBiμ,detBκ0

nur aus der Null besteht. Für den Falle eines Polynoms der Form B1μ+B0 ergibt sich daher die folgende Überlegung. Man beachte, daß man auf dieses lineare Polynom alle anderen Polynome zurückführen kann, u.U. durch Vergrößern der Dimension, z.B. durch das Begleitpolynom.

Im monischen Falle, also detB10, besteht das Spektrum des besagten linearen Polynomes B1μ+B0 genau dann nur aus der Null, wenn B11B0 ähnlich ist zu

diag(0,,0,J1,,Jk),Ji=(δp+1,q)p,q=1sν,

also die Matrizen Ji Jordankästen zum Eigenwert Null sind. Der Fall detBκ0 ist stets zu erreichen, wenn alle Stufen implizit, oder sogar blockimplizit sind.

Bei expliziten Stufen oder bei Mischung aus expliziten und impliziten Stufen, ist das reguläre, oder häufig sogar singuläre Matrizenbüschel B1μ+B0 entsprechend zu untersuchen, siehe Gantmacher, Felix R. (1908--1964), Gantmacher (1986).

_{Tischer, Peter E.}Tischer (1983)2 und {Tischer, Peter E.}{Sacks-Davis, Ron}Tischer/Sacks-Davis (1983)3 z.B. gehen hier den anderen möglichen Weg mit B00.

Dennoch bedeutet B0=0 im Falle B1μ+B0, bzw. Bi=0 (i=0,,κ) im Falle i=0κBiμi, einen großen Rechen- und Speichervorteil, wenn der Prädiktor entsprechend gewählt wird. Sehr vereinfacht ausgedrückt heißt dies: ``Je mehr Nullen in den Stufen, desto besser."\space Die Vereinfachung liegt hier maßgeblich u.a. daran, daß der Prädiktor entsprechend gewählt werden muß und, daß Nullen bei den ym+i-Werten nicht automatisch auch zu Nullen bei der Implementierung führen, welche ja auf der Darstellung in Form von rückwärtsgenommenen Differenzen beruht. Man beachte, daß bei z.B. 1000 Schritten u.U. 2000--9000 mehr Vektoradditionen und Skalar-Vektor-Multiplikationen ausgeführt werden müssen, wenn man “Nullen verschenkt.”\space Es ist nicht ungewöhnlich, daß eine Funktionsauswertung billiger ist als die Bildung einer längeren Linearkombination von gespeicherten Vektoren. % _{Adams, J.C.}% Adams-Formeln und BDF sind hier typische Vertreter für beide möglichen Extrema, was z.T. ihren hohen Grad an Effizienz plausibel macht.

4. Beispiel: Für das explizite Eulerverfahren yn+1yn=hfn erhält man mit den Vektoren

un+1=(ynyn+1)un=(yn1yn)

die vier Matrizen

A0=(0100),A1=(1011),B0=(0100),B1=(1010).

Die Matrizen A1, A0, B1 und B0 sind natürlich nicht eindeutig. An der Nichtinvertierbarkeit von B1 ändert dies aber nichts. Offensichtlich hat das Matrixpolynom B1μ+B0, also hier das verallgemeinerte Eigenwertproblem, einen Eigenwert bei unendlich und damit ist das Verfahren natürlich nicht A0-stabil. An dieser Stelle werden auch die Überlegung bezüglich des singulären Matrizenbüschels von Wichtigkeit.

5. Es sei nocheinmal darauf hingewiesen, daß für die Ordnung p=1 und p=2 keine originär neuen Formeln gesucht und damit auch nicht gefunden werden konnte. Es wurde tatsächlich dreimal hintereinander dasselbe Verfahren bei den beiden niedrigsten Ordnungen als zyklisches Verfahren benutzt. Im Hinblick darauf, daß sich die beiden Stabilitätsmerkmale, wie _{Widlund, Olof B.}

Widlund-α-Winkel, also A[α]-Stabilität und Widlund-Distanz, also S[δ]-Stabilität, nicht weiter verbessern lassen, wird dies verständlich. Zudem bleibt damit auch der Wert γ=\frac{β}/{α}, der bei der Iterationsmatrix W=IhγJ wichtig wird, über den gesamten Zyklus konstant, obwohl dies nicht so wesentlich ist, wie vermutet werden könnte. Kleinere Fehlerkonstanten und ^{Schrittweitenwechsel-Stabilitätseigenschaften} blieben bei der Auffindung der neuen zyklischen Formeln unberücksichtigt.

6. Das Tschebyscheffsche-α-Äquilibrierungsmaß μ einer linearen Mehrschrittformel ist definiert durch

μ=1|ακ|maxi=0κ(|αi|,|βi|).

Hiermit wird ein erster Eindruck von der unterschiedlichen Gewichtung einzelner Terme in dem Verfahrensausdruck vermittelt.

Die zyklischen Formeln von Tendler weisen annähernd ein Äquilibrierungsmaß von 1:p auf, wobei p die Ordnung des Zykluses ist.

Die BDF hingegen weisen annähernd ein Äquilibrierungsmaß von 1:p2 auf, welches also lediglich halb so groß ist, wie die Äquilibrierungsmaße der Formeln von J.M. Tendler.

Diesem Punkte wurde beim Aufsuchen der Formeln keinerlei Berücksichtigung zuteil. Anders ist dies etwa bei den Formeln von Filippi/Kraska (1973). {Filippi, S.}{Kraska, E.}

Hier wurden die Formeln durch Lösungen eines linearen Optimierungsproblems bestimmt, wobei die Summe der Beträge der Koeffizienten der Verfahren in der Zielfunktion auftauchten. Die Formeln von {Tischer, Peter E.}Tischer (1983)2 und Tischer/Sacks-Davis (1983)3{Tischer, Peter E.}_{Sacks-Davis, Ron} weisen besonders ungünstige Äquilibrierungsmaße auf. Die Maße bewegen sich von ca. 1:p2 bis zu ca. 1:2p2.

7. Zusammenstellend noch die Werte

γ=maxi=1βiiαii,

welche wichtig sind für die Konvergenz der Picard-Iteration. Die Werte sind auf eine Dezimale hinter dem Komma genau.

p1234567γ1.00.70.50.50.50.50.4

8. Es seien jetzt wieder nur die Tendlerschen Zyklen betrachtet. Die Matrix Q sei diejenige Matrix, wie sie bei der Umwandlung von Ordinatenwerten in Werte für rückwärtsgenommene Differenzen auftaucht. Durch Vormultiplikation mit der Matrix VQV, passenden Ranges, mit V=(δn+1i,j), ergeben sich für die Ordinatenwerte ym+i die Werte in der folgenden Tabelle, also die Werte für die Darstellung der Formeln in Form von rückwärtsgenommenen Differenzen. Dieses sind jedoch noch nicht die Werte, die nachher wirklich gespeichert werden. Die Werte zm+i bleiben hier unberücksichtigt, da sie nicht in der Form rückwärtsgenommener Differenzen verwendet werden. Dies liegt an der günstigen Eigenschaft der Tendlerschen Zyklen, daß die skalierten Ableitungswerte niemals interpoliert zu werden brauchen. Bei einem Schrittweitenwechsel muß ein Ableitungswert lediglich neu skaliert werden. Dies ist eine typische Eigenschaft der Tendlerschen Zyklen. Für anders gebaute Zyklen gilt dies nämlich nicht so ohne weiteres.

pi23456711,2,3121,2,32131,2632369041,21264336013841151,26030201512354039018014511842760376011011513361603020151210211406304103052432023720126027027022919547380861015030532928571,2420210140105847060312602430510405313252210454607350364013651148931774

Wegen der ersten Zeile der Matrix Q und der Konsistenzbedingung ρ(1)=0, wird der Wert für 0 nicht mit vermerkt. Die erste Stufe bei allen zyklischen Formeln ist die BDF und hier gilt

i=1κ1iiyn+1=zn+1.

Die Normierung von yn+1 wird hier nicht gesondert durchgeführt.

3. Die Fehlerfaktoren der Korrektorformeln

The new methods presented were obtained solely on the basis of optimizing the Widlund wedge angle associated with the method. It is felt that attention should also be focused on simultaneously minimizing the local discretization error. In addition, higher stiffly stable methods should be derived.

J.M. Tendler (1973)\silentreference{Tendler (1973)}

Nun seien die Beträge der Fehlerfaktoren der sieben zyklischen Verfahren angegeben. Zum Vergleich werden auch die Fehlerfaktoren der BDF angeschrieben. Diese sind gerade die Fehlerfaktoren der ersten Stufe der 7 zyklischen Verfahren. Mit Fehlerfaktor ci,p+1 wird hier die erste nichtverschwindende Komponente des Matrix-Vektor-Produktes aus Cp+1,k und den Koeffizienten (α,β) des Verfahrens bezeichnet, skaliert mit einer entsprechenden Fakultät. Es ist also

ci,p+1:=1(p+1)!Cp+1,k(αβ).

1. Die Fehlervektoren γp für alle Ordnungen p=1,,7 und alle Stufen, lauten exakt ausgerechnet wie in der Tabelle sichtbar.

p1234567c1,p+1122332125106071052c2,p+112233212510121071052c3,p+1122312109911827251514c4,p+123921699713192

2. Normiert man stufenweise den führenden Koeffizienten des Verfahrens auf eins, so normieren sich entsprechend die Werte ci,p+1, für i=1,,. Für die zyklischen Formeln von Tendler ergeben sich dann nach kurzer Rechnung, die Werte in der entsprechenden Tabelle.

p1234567c1,p+1/α110.50.20.1360.0960.072990.058310.04821c2,p+1/α220.50.20.1360.0960.072990.059000.04821c3,p+1/α330.50.20.160.10750.072100.057360.06281c4,p+1/α440.10320.10560.067654

3. _{Henrici, Peter Karl Eugen} Zur Gegenüberstellung nun die gleichen Fehlerfaktoren diesmal jedoch skaliert mit dem Kehrwert der Summe der β-Koeffizienten, was gerade zur Fehlerkonstante von Henrici führt. Man beachte, daß hier erneut stets stufenweise die Fehlerfaktoren gebildet wurden.

p1234567c1,p+1/σ1(1)0.50.30.250.20.160.142860.125c2,p+1/σ2(1)0.50.30.250.20.160.151630.125c3,p+1/σ3(1)0.50.30.0830.160.1830.234520.14257c4,p+1/σ4(1)0.04330.032890.12079

Auffällig ist, daß die Tabellenwerte betragsmässig leicht äquilibrierter erscheinen und im ganzen gesehen, etwas größer geworden sind. Natürlich gilt hier immer σip(1)0, für alle Stufen i=1,, und alle Ordnungen p=1,,7. Das heißt, 1 ist einfache Wurzel jeder Stufe, nicht notwendig die betragsmässig größte.

4. Übersichtlich zum Vergleich gegenübergestellt, werden jetzt die zyklusweise gebildeten Fehlerkonstanten von Henrici aufgeführt. Man beachte, daß diese Konstante jeweils immer zyklusweise gebildet wird. Die Werte lauten für alle Ordnungen p=1,,7:

p1234567Cp1.513.751.41914891.6872031.63206970.684573

Diese Konstanten sind nun betragsmässig am größten, von allen bisher angegebenen Konstanten. Auffallend ist hier die Konstante für die Ordnung p=3.

5. Um einen Überblick über die Größenordnungen zu erhalten, seien auch die ^{Fehlerkonstanten der Adams-Moulton-Verfahren} {Adams, J.C.}{Moulton, F.R.} in der folgenden Tabelle angeschrieben.

p1234567Cp0.50.8330.04170.02640.01880.01430.0114

Wegen ρ(1)=σ(1)=1 und αk=1 stimmen alle drei verschiedenen Arten an Fehlerkonstanten über-ein. Vergleicht man nun alle bisher angegebenen Konstanten, so muß man beachten vor welcher h-Potenz diese Konstanten stehen. Beispielsweise führt bei einem Verfahren siebenter Ordnung eine Halbierung der Schrittweite asymptotisch betrachtet, zu einer Erhöhung der “Genauigkeit” um den Faktor 256. Überdies ist der Unterschied zwischen einem Verfahren p-ter und (p+1)-ter Ordnung nur gering, falls die Schrittweite variabel sein kann. Ist das Verfahren p-ter Ordnung nicht mehr ausreichend genau genug, so ist es gerade die zentrale Aufgabe der Ordnungssteuerung dies zu erkennen und durch Wahl einer passenderen Ordnung, den Genauigkeitsverlust wieder auszugleichen. Weiterhin beachte man die Wachstumseigenschaften der Fehlerkonstante von Henrici bei mehrfacher Wiederholung einer Stufe.

4. Die Stabilitätseigenschaften der Korrektorformeln

While methods suitable for non-stiff problems, , are easy to find and to analyse, methods appropriate for stiff problems seem to be able to hide their presence much more successfully. However, there are such methods

J.C. Butcher (1987) {Butcher, John Charles}

Zum Vergleich seien die Werte für den Widlund-α-Winkel der Formeln von J.M. Tendler mit denen der BDF in der untenstehenden Tabelle gegenübergestellt. Man beachte, daß die BDF7 nicht mehr D-stabil ist, erst recht also nicht mehr A[α]-stabil. Rechts daneben ist auch eine Gegenüberstellung der Widlund Distanzen δ, bei der S[δ]-Stabilität.

Widlung-α-Winkel.

p TENDLER BDF
1 90.00 90.00
2 90.00 90.00
3 89.43 86.03
4 80.88 73.35
5 77.48 51.84
6 63.25 17.84
7 33.53 *

Widlung-Distanz δ.

p TENDLER BDF
1 0.0000 0.0
2 0.0000 0.0
3 0.0048 0.083
4 0.2400 0.67
5 1.4000 2.3
6 2.9000 6.1
7 10.2000 *

Die Definition der Begriffe Widlund-Winkel, Widlund-Distanz, A[α]-Stabilität, S[δ]-Stabilität findet man bei Tendler (1973) und Tendler/Bickart/Picel (1978), {Tendler, Joel Marvin}{Bickart, Theodore A.}{Picel, Zdenek}% aber sinngemäß auch bei Albrecht (1979),{Albrecht, Peter} §7.2., und weiteren Autoren. \else Betrachtet man die BDF wie üblich als einstufige Verfahren, so ergeben sich leicht größere Widlund-Winkel für die BDF\null. _{N\o rsett, Syvert Paul}N\o rsett (1969) stellte die folgende Tabelle der Widlund-α-Winkel für die BDF auf, hier als einstufige Verfahren betrachtet und gab an, daß der Fehler unter 1/60=0.016 liege.

Ordnung 1 2 3 4 5 6
αmax 90 90 88.45 73.23 51.83 18.783

5. Die sieben Prädiktorformeln

\beginepigram_{Butcher, John Charles}

The difficulty of resolving various other questions concerned with incorporating general linear methods into practical software should not be underestimated.

\author J.C. Butcher (1987) \endepigram

Es seien jetzt noch die Prädiktorformeln angegeben. Bei ihnen handelt es sich um nichts anderes als die expliziten BDF\null. Diese Formeln sind für Ordnungen p>2 nicht mehr D-stabil, d.h. sie dürfen nicht alleinstehend ohne einen Korrektor verwendet werden. Das Verfahren erster Ordnung ist natürlich nichts anderes als das explizite Eulerverfahren yn+1=yn+hfn, und die Prädiktorformel zweiter Ordnung ist die (explizite) Mittelpunktsregel yn+1=yn1+2hfn (= explizite ^{Nyström-Formel zweiter Ordnung}).

Die Angabe der Werte geschieht zunächst in der Form

yn+10=αnyn++αnp+1ynp+1+βnzn.

für die Ordnungen p=1,,7. Die Werte βni für i>0 sind sämtlich Null, d.h. heißt es werden vom Prädiktor keine weiteren vergangenen skalierten Ableitungen benutzt.

Zuerst nun die Werte αni und βn in der obigen Darstellung mit gewöhnlichen Ordinatenwerten ym+i, siehe erste Tabelle.

kznynyn1yn2yn3yn4yn5yn611122013332621244103183631355651212012601220123126677101501010010501015102107760960126060105060700603156084601060

Die Darstellung der Prädiktorformeln in Form von rückwärtsgenommenen Differenzen, also in der Darstellungsform

yn+1=α~n0yn++α~np+1p1yn+βnzn

für die Ordnungen p=1,,7, entnimmt man der zweiten Tabelle.

kzn023456111221133121244132213551432231466154233241577165243342516

In kompakter Form seien noch diejenigen Werte angegeben, die auch tatsächlich im Programm TENDLER benutzt werden, wenn auch nicht in dieser kompakten Form. An späterer Stelle wird nocheinmal kurz auf die Herleitung dieser Koeffizienten eingegangen. Die nachstehende Tabelle enthält in gewisser Hinsicht ein “Gemisch” von Prädiktor- und Korrektorformeln. Dabei werden für y rückwärtsgenommene Differenzen verwendet, für z jedoch nicht.

ki23456z1z0z111,2,3121,2,32331,292541123315234132241,222383171825343994783545451,212512101247136374813712532532410012407073601231601373120110545742581710169240615311012061,2137101172046151911201973001471063353255711001819600795039196000147710072064352920018894001609600352724009311500307920017532071,234320303402536058924010175234036320732661522130749180733080360010318054730127413167511511503689900121504003300025745027371501215

Hierbei beachte man, daß die Werte in Tendler/Bickart/Picel (1976) {Tendler, Joel Marvin}{Bickart, Theodore A.}_{Picel, Zdenek} für p=6, i=4 nicht sämtlich richtig sind. In dem Programm STINT und auch in dem Programm DIFJMT werden jedoch die korrekten Werte gespeichert.