Zur Lösung des Differentialgleichungssystems 12.13 werden in JANAP die ``Backward-Differentiation Formeln'' (BDF) mit automatischer Schrittweiten- und Ordnungssteuerung verwendet [79][70][69][68][67][66][65][37][17].
Das Prinzip der BDF-Methode besteht aus dem folgenden Algorithmus:
Wir gehen von unserer Problemstellung 12.13 aus:
Die zentrale Aufgabe ist nun die Ermittlung eines Wertes für
für den nächsten Zeitschritt,
also die Lösung von
Als Basis liegt uns das Ergebnis der
Integration unseres Differentialgleichungssystems 12.13
für die letzten 
Zeitschritte vor, also für die Zeitpunkte
.
kann gemäß den BDF-Formeln für die -te Ordnung
mit Formel 12.16 approximiert werden.
wobei die 
Konstante sind und 
die Größe des nächsten
Zeitschritts (
) ist.
Wenn wir die approximierten Werte in 12.15 einsetzten, erhalten
wir das nichtlineare algebraische Gleichungssystem 12.17.
Gleichungssystem 12.17 kann nun mit der in 12.1 dargestellten Methode gelöst werden.
Dieser Vorgang wird solange wiederholt, bis der angegebene Endzeitpunkt der transienten Analyse erreicht ist.
Die Schrittweite ist variabel, um das steife Differentialgleichungssystem
effizient zu lösen.
Um 
nach Gleichung 12.16 bestimmen
zu können, müssen die Werte 
in dieser Gleichung
vorgegeben sein.
Da die gleichen 
für jedes Element von
verwendet werden, reicht es, wenn wir 
für ein
Element des Vektors ableiten.
Wir haben daher die 
für die folgende Beziehung zu ermitteln:
Die Strategie ist, daß die 
derart ermittelt werden,
daß Gleichung 12.18 das exakte Ergebnis liefert, wenn
ein Polynom vom Grad ist, also
gilt.
Wir können die Werte von der letzten 
Zeitschritte
in Form eines linearen Gleichungssystems mit Rang zusammenstellen:
Gleichung 12.19 kann auch als Summe von Polynomen der
Ordnung dargestellt werden, wobei das Polynom 
genau zum Zeitpunkt 
den Wert 1 hat, zu den anderen Zeitpunkten
jedoch 0 ist.
Damit gilt
Durch das Einsetzen von 
in Gleichung 12.22 ergibt
sich
Durch den Vergleich der Koeffizienten von 12.18 und 12.23 ergibt sich
Die 
können durch das Lösen des linearen Gleichungssystems
12.25 berechnet werden.
Der Rang des Gleichungssystems beträgt .
Gleichungssystem 12.25 kann a priori gelöst werden wodurch
folgende Formeln für 
entstehen:
Die 
müssen für jeden Zeitschritt und bei jeder Änderung
der Ordnung neu berechnet werden.
Da dies sehr zeitaufwendig ist, wurden abgekürzte Formeln entwickelt [17],
die die 
aus den Werten des letzten Zeitschritts bzw.
der alten Ordnung errechnen.
Ein weiteres Problem stellt die Genauigkeit der Berechnung dar.
Die Formeln 12.26 verwenden die letzten Zeitpunkte
zur Berechnung. Des weiteren verwendet Formel 12.16
die 
-Werte der letzten Zeitpunkte. Diese Werte sind in Relation
zu den Änderungen der -Werte sehr groß wodurch Instabilitäten
bei der numerischen Auswertung von 12.16 während der
Simulation auftreten können.
Eine numerisch stabilere Beziehung zur Berechnung von
ist
wobei
gilt.
Sowohl als Startwert für die Newton-Iteration (12.1) als
auch zur Abschätzung des Fehlers und die darauf aufbauende
Schrittweiten- und Ordnungssteuerung wird ein Schätzwert für
benötigt.
Dieser Wert kann nicht einfach durch Integration von 12.16
ermittelt werden, da bei der Integration ein Polynom von Ordnung
entstehen würde. Um alle Koeffizienten dieses Polynoms zu
berechnen, wären dann 
Werte notwendig - wir haben aber nur
Werte. Daher wird wie bei der Berechnung der Ableitung
ein Polynom mit Ordnung verwendet:
Ähnlich wie die 
sind die 
die Lösung eines
linearen Gleichungssystems 12.30 mit Rang .
Die kompakte Darstellung der lautet:
Auch hier gibt es wieder eine numerisch stabilere Variante, die auf den Differenzen der x-Werte beruht:
wobei
gilt.
In praktisch allen in diesem Kapitel abgeleiteten Gleichungen kommt
die Schrittweite 
als zentrales bestimmendes Element
vor. Wie groß ist bzw. soll die Schrittweite gewählt werden.
Die oberste Regel lautet:
Die Schrittweite soll so groß wie möglich gewählt werden, sodaß die geforderte Genauigkeit noch erreicht wird und die Konvergenz des Verfahrens gewährleistet werden kann.Aus der Wahl der Berechnung von
als Polynom vom
Grad ergibt sich, daß das Integrationsverfahren nur genau dann
exakt ist, wenn 
selber ein Polynom mit einem Grad
von oder kleiner ist.
Sofern diese Bedingung nicht erfüllt ist - das ist der Regelfall -
muß der Fehler durch die Diskretisierung der Integration des
Differentialgleichungssystems abgeschätzt werden und mit der
vorgegebenen Genauigkeit verglichen werden.
Wenn diese Genauigkeit nicht erreicht wurde, so kann man davon ausgehen,
daß die Schrittweite zu groß gewählt worden ist.
Man muß daher diesen Schritt mit einer kleineren Schrittweite wiederholen.
Ist die Genauigkeit erreicht oder wesentlich besser als gefordert, kann
man die Schrittweite (für den nächsten Schritt)
entsprechend erhöhen.
Der ganze Mechanismus steht und fällt natürlich mit der Bestimmung
des Fehlers.
Ideal wäre es, den exakten Fehler ermitteln zu können. Dazu muß jedoch
die exakte Lösung bekannt sein - die wollen wir aber ausrechnen.
Statt dessen müssen wir eine Abschätzung des Fehlers verwenden.
Ein mit vernünftigen Aufwand zu berechnender Fehler ist der
``lokale Diskretisierungsfehler'' (local truncation error),
das ist der Fehler, der durch die Verwendung eines Polynoms -ter
Ordnung bei einem Zeitschritt entsteht.
Der lokale Diskretiersierungsfehler ist definiert zu
also die Abweichung des errechneten Wertes 
vom
Sollwert 
. Wie schon vorher erwähnt, ist aber
der Sollwert nicht bekannt.
Wir schätzen den Sollwert daher mit dem vor Lösung des
Differentialgleichungssystems ermittelten Wert 
nach
Gleichung 12.29 oder 12.32 ab.
Der Fehler dieser Abschätzung muß in der Größenordnung von
liegen, wenn wir zur Approximation ein Polynom der Ordnung
verwendet haben.
Damit gilt
nach Gleichung 12.36 wird als Abschätzung des
lokalen Diskretisierungsfehlers herangezogen.
Bei der konkreten Implementierung ist zu beachten, daß das ansich
stabile BDF-Verfahren instabil werden kann, wenn sich die Schrittweite
rasch ändert. Die Änderung der Schrittweite muß daher beschränkt werden.
Ein weiterer wichtiger Parameter während der Integrationsphase ist
die gewählte Ordnung des zur Approximation von 
verwendeten Polynoms.
Es werden folgende Kriterien zur Wahl der Ordnung verwendet:
Koeffizienten des Polynoms 12.16 vorhanden sind.
Wie in [119][70] gezeigt wurde, können algebraische Differentialgleichungssysteme im Gegensatz zu den gewöhnlichen Differentialgleichungssystemen im allgemeinen nicht mit den BDF-Formeln gelöst werden, wenn automatische Schrittweiten- und Ordnungssteuerung verwendet wird.
Bei den Netzwerkgleichungen handelt es sich jedoch immer um gewöhnliche Differentialgleichungen, sofern der Differentialoperator nur in den Formen
vorkommt. Diese Bedingung ist bei allen idealen Bauelementen und auch bei der Modellierung von Halbleiterbauelementen [100] erfüllt. Falls Simulationen durchgeführt werden sollen, die diese Bedingung nicht einhalten, sind geeignete Maßnahmen zu treffen (Beschränkung der Ordnung, Fixierung der Schrittweite) [70].
Eine besonders unangenehme Situation ist das Auftreten von Unstetigkeitsstellen
in der Lösung des Differentialgleichungssystems.
Diese Unstetigkeitsstellen werden durch die in JANAP verfügbaren
Bauelemente (Schalter) oder Methoden zur Modellierung
(einige eingebaute Funktionen wie IFEQ) zusätzlich hervorgerufen.
Die Reaktion des Integrationsverfahrens bei einer Unstetigkeitsstelle
ist die folgende:
In [69] wird eine Methode aufgezeigt, bei der eine Unstetigkeitsstelle schneller (in weniger Schritten und mit einer größeren Ordnung) verarbeitet werden kann.
