Differenzengleichung (Differenzenverfahren)

Grundidee und Begriffe

Eine Differenzengleichung ist eine rekursive, numerisch lösbare Berechnungsvorschrift für eine diskret definierte Folge nummerierter Werte. Diese Werte heißen Folgeelemente, Stützstellen oder Knoten und werden meist als y_k bezeichnet. Zwischen zwei benachbarten Werten liegt ein konstantes Intervall Δx oder bei zeitabhängigen Systemen eine Schrittweite Δt = h.

Differenzengleichungen werden in vielen Bereichen verwendet, etwa Wirtschaft, Medizin, Technik, Elektrotechnik, Regelungstechnik, Kybernetik, Informatik und Akustik. Sie sind besonders wichtig, weil viele Differenzialgleichungen nicht einfach geschlossen lösbar sind. Dann ersetzt man den Differenzialquotienten durch einen Differenzenquotienten. Dadurch entsteht automatisch eine rekursive Gleichung: Ein aktuelles Folgeglied hängt von einem oder mehreren vorherigen Folgegliedern ab.

Eine Differenzengleichung heißt auch Rekursionsgleichung. Sie ist von n-ter Ordnung, wenn die höchste Ordnung der vorkommenden Differenzen gleich n ist. Bei einer Differenzengleichung 1. Ordnung bezieht sich jedes y_k auf ein zurückliegendes y_{k-1}. Bei höherer Ordnung können mehrere frühere Werte beteiligt sein, zum Beispiel y_k = f(x_k, y_{k-1}, y_{k-2}, y_{k-3}, ...). Die Indizes laufen typischerweise über k = [0,1,2,3,...,k_max].

Die numerische Lösung ergibt meist eine Tabelle. Für zeitabhängige Systeme enthält sie etwa k, die Zeit k·Δt, die Ausgangsfolge y_k und eine Eingangsfolge u_k. Die Eingangsfolge u_k ist nicht rekursiv; sie kann zum Beispiel eine normierte Sprungfunktion u_k = 1(t_k), die Ausgangsfolge eines vorgeschalteten Systems oder eine Tabelle sein. Die Anzahl der Folgen ergibt sich aus Anzahl der Folgen = t_max / h + 1. Dabei ist t_max der Darstellungszeitraum und h die Schrittzeit. Je kleiner h ist, desto genauer wird die Approximation, aber desto mehr Folgeglieder müssen berechnet werden.

Einfache Folgen und Anfangswertprobleme

Bei einfachen Aufgaben wie Zinseszins, Bevölkerungswachstum oder dem Füllen und Entleeren von Behältern liegt nicht unbedingt zuerst eine Differenzialgleichung vor. Man startet mit einem Anfangswert y_0 und berechnet für Zeit t, Anzahl der Folgeglieder k_max und Schrittweite h die Folgewerte y_k. Aus jedem früheren Folgeglied wird das nächste berechnet.

Wichtige einfache Typen sind arithmetische und exponentielle Folgen. Bei einer arithmetischen Folge wächst oder fällt jedes Folgeglied um einen festen Betrag. Bei einer exponentiellen Folge wächst oder fällt jedes Folgeglied um einen relativen Anteil. Eine typische Gleichung lautet y_{k+1}=y_k ± C·y_k·h = (1 ± C·h)·y_k. Für eine beliebige Folge k kann man direkt schreiben: y_k = (1 ± C·h)^k · y_0. Diese Form entspricht auch der Zinseszinsformel K = (1 + p/100 · 1 Jahr)^{n Jahre} · K_0. Der Wachstumsfaktor ist q = K_{k+1}/K_k = (1 + p/100).

Beim Anfangswertproblem sucht man die spezielle Lösung einer gewöhnlichen Differenzialgleichung unter Berücksichtigung eines Anfangswerts. Für eine Gleichung y'(x)=f(x,y) wird y(x_0)=y_0 vorgegeben, meist für x=0. Wird die Ableitung durch den Vorwärts-Differenzenquotienten ersetzt, gilt näherungsweise y'(x) ≈ (y_{k+1}-y_k)/h. Daraus entsteht die explizite Differenzengleichung (y_{k+1}-y_k)/h = f(x_k,y_k) und damit y_{k+1}=y_k+f(x_k,y_k)·h. Dies entspricht dem Euler-Vorwärts-Verfahren. Kleinere Schrittweiten h verringern die Integrationsfehler. Verfahren wie das Trapezverfahren von Heun oder das Runge-Kutta-Verfahren erlauben bei gleicher Genauigkeit größere Schrittweiten.

Differenzenquotienten und dynamische Systeme

Gewöhnliche lineare Differenzialgleichungen dynamischer Systeme können durch Austausch der Differenzialquotienten gegen Differenzenquotienten in Differenzengleichungen überführt werden. Ein Beispiel für ein System 1. Ordnung ist a_1·y'(t)+a_0·y(t)=b_0·u(t). Die Eingangsgröße ist u(t), die Ausgangsgröße y(t), die unabhängige Variable ist die Zeit t.

Der Vorwärts-Differenzenquotient lautet y'(t) ≈ (y_{k+1}-y_k)/h. Für differenzierende Systeme wird entsprechend u'(t) ≈ (u_{k+1}-u_k)/h verwendet. Der Rückwärts-Differenzenquotient lautet y'(t) ≈ (y_k-y_{k-1})/h und für das Eingangssignal u'(t) ≈ (u_k-u_{k-1})/h. Der zentrale Differenzenquotient lautet y'(t) ≈ (y_{k+1}-y_{k-1})/(2h). Beim zentralen Differenzenquotienten steigt die Genauigkeit nicht nur mit fallendem h, sondern mit dem Quadrat des fallenden h.

Lineare dynamische Systeme werden oft als Übertragungsfunktion G(s) beschrieben. Für phasenminimale Elementarsysteme ergeben sich durch inverse Laplace-Transformation Differenzialgleichungen 1. und 2. Ordnung. Genannt werden unter anderem I-Glied, D-Glied, PD_1-Glied, PT_1-Glied, PT_{2kk}-Glied und Totzeitglied. Ein Totzeitglied y(t)=K·u(t-T_t) entsteht durch Laufzeiten von Material oder Signalen; seine Laplace-Transformierte ist keine gebrochene rationale Funktion, lässt sich numerisch aber leicht behandeln.

Bei der Herleitung für lineare Systeme unterscheiden sich Vorwärts- und Rückwärtsverfahren vor allem in der zeitlichen Zuordnung der Werte. Beim Vorwärtsverfahren wird mit y_0 gestartet, und y_1 sowie alle weiteren Werte werden rekursiv berechnet. Beim Rückwärtsverfahren wird y_k aus y_{k-1} berechnet; es kann bei k=0 bereits ein kleiner Wert entstehen, obwohl der Anfangswert eines Systems eigentlich Null sein sollte. Für ein PT_1-Glied ergibt sich mit Vorwärtsquotient y_{k+1}=y_k+(K_{PT1}·u_k-y_k)·Δt/T_1. Mit Rückwärtsquotient lautet eine vereinfachte Form y_k=y_{k-1}+(K_{PT1}·u_k-y_{k-1})·Δt/(T_1+Δt).

Beispiele aus Wachstum, Bewegung und Schwingung

Beim Bevölkerungswachstum wird ein exponentielles Wachstum berechnet. Gegeben sind B_0 = 20 Millionen, ein konstantes Bevölkerungswachstum G = 4 % pro Jahr, h = 1 Jahr und der Anfangswert B_0 = 20 Millionen. Die Differenzengleichung lautet B_{k+1}=B_k+B_k·a·h = B_k + B_k·0,04 = B_k·1,04. Der Faktor 1,04 ist der Wachstumsfaktor pro Folgeglied. Nach 50 Jahren ergibt sich B_{50}=20·1,04^{50}=142,1336669 Millionen. Allgemein gilt B_k=B_0·(1+G/100)^k. Der Artikel weist darauf hin, dass ein solches ungebremstes Wachstum in der Praxis nicht dauerhaft vorkommt, weil andere Einflüsse wie Nahrungsmittelknappheit dagegen wirken.

Ein Beispiel für ein Anfangswertproblem ist y'(x)=y(x)+e^x mit analytischer Lösung y(x)=(x+y_0)·e^x. Für y_0=5 und h=0,05 entsteht durch den Vorwärts-Differenzenquotienten (y_{k+1}-y_k)/h = y_k + e^{x_k}, also y_{k+1}=y_k+(y_k+e^{x_k})·h. In der Tabelle werden für n=50 Folgen numerische Werte mit der analytischen Lösung verglichen; bei x=2,5 steht für das Differenzenverfahren y_{50}=85,446, für die analytische Funktion y_{50}=91,369.

Beim Fallschirmspringer lautet die nichtlineare Bewegungsgleichung mit Luftwiderstand m·v̇(t)=m·g-c·v^2(t). Mit c=0,32 kg/m, m=80 kg, g=10 m/s² und Anfangswert v_0=0 erhält man v_{k+1}=v_k+(g-c/m·v_k^2)·h. Für h=1 s lauten die ersten Werte v_1=10, v_2=19,6 und v_3=28,063 m/s. Die Folge nähert sich asymptotisch 50 m/s; ab v(t)=50 m/s nimmt die Fallgeschwindigkeit nicht mehr zu. Mit h=0,01 s wird nach 13 s eine Fallstrecke von etwa 478 m angegeben.

Beim mathematischen Pendel lautet die Bewegungsgleichung φ̈ + (g/l)·sin(φ)=0. Durch ω=φ̇ wird daraus ein System erster Ordnung: ω̇=-(g/l)·sin(φ) und φ̇=ω. Numerisch wird ω_{k+1}=ω_k-(g/l)·sin(φ_k)·h berechnet. Für den Winkel wird φ_{k+1}=φ_k+ω_{k+1}·h verwendet. Die Nutzung der zuvor berechneten Winkelgeschwindigkeit verbessert die Genauigkeit; eine nochmalige Anwendung des expliziten Euler-Verfahrens würde eine aufklingende, instabile Lösung ergeben.

Höhere Ordnung, Zustandsraum und Simulation

Für dynamische Systeme mit Differenzialgleichungen zweiter und höherer Ordnung nennt der Artikel drei Wege: Lösung einer DGL 2. Ordnung mit einem Modellregelkreis, Lösung einer DGL höherer Ordnung mit Differenzenquotienten und Lösung mit Zustandsvariablen nach der Zustandsraumdarstellung. Bei einem Schwingungsglied PT2_{kk} wird zum Beispiel T^2 y''(t)+2DTy'(t)+y(t)=K·u(t) verwendet.

Bei der direkten Methode mit Differenzenquotienten braucht jede Ableitung einen passenden Differenzenquotienten. Für die zweite Ableitung lauten die Formen d²y/h² ≈ (y_{k+2}-2y_{k+1}+y_k)/h² im Vorwärtsverfahren und d²y/h² ≈ (y_{k-2}-2y_{k-1}+y_k)/h² im Rückwärtsverfahren. Im Beispiel G(s)=1/(0,06s²+0,1s+1) gehört die Differenzialgleichung 0,06·y''(t)+0,1·y'(t)+y(t)=u(t). Mit h=0,01 und Sprung u_k=1 ergibt sich bei der dargestellten Rechnung der Maximalwert der ersten Amplitude bei k=78 und t=0,78 s mit y(t)=1,5192.

In der Zustandsraumdarstellung wird eine Differenzialgleichung n-ter Ordnung in n gekoppelte Zustands-Differenzialgleichungen erster Ordnung überführt. Zustandsvariablen beschreiben physikalisch den Energiegehalt der Speicherelemente eines dynamischen Systems. Für das Beispiel 0,06·y''(t)+0,1·y'(t)+y(t)=u(t) wird y''(t)=16,666·u(t)-1,6666·y'(t)-16,666·y(t). Numerisch verwendet der Artikel unter anderem y''k=16,666·u_k-1,666·x2_k-16,666·x1_k, x2{k+1}=x2_k+y''k·h und x1{k+1}=x1_k+x2_{k+1}·h. Mit h=0,05 wird der Maximalwert der ersten Amplitude bei k=15 und t=0,75 s mit ungefähr y(t)=1,515 angegeben. Gegenüber der Methode mit Differenzenquotienten 2. Ordnung gilt die Regelungsnormalform bei gleicher Genauigkeit als deutlich überlegen, weil weniger Folgeglieder nötig sind.

Differenzengleichungssysteme erlauben auch die Simulation mehrerer zusammenwirkender Systeme, zum Beispiel Reihenschaltungen, Parallelschaltungen oder Rückführungen. Ausgänge eines Systems werden dabei Eingänge des nächsten. Lineare Verzögerungsglieder wie PT1 können durch ideale differenzierende Glieder wie PD1 bei gleichen Parametern kompensiert werden. Im Beispiel mit T=1, K=1 und h=0,2 wird ein zweites PT1-Glied durch ein ideales PD1-Glied vollständig kompensiert.

Nutzen, Grenzen und Sonderfälle

Ein Vorteil des numerischen Differenzenverfahrens ist, dass es gegenüber der konventionellen Lösung von Differenzialgleichungen relativ geringe mathematische Kenntnisse erfordert. Differenzengleichungen können direkt über Differenzenquotienten aus Differenzialgleichungen gebildet werden und lassen sich mit jeder Programmiersprache oder auch mit Tabellenkalkulation berechnen.

Der Approximationsfehler beim Vorwärts- oder Rückwärtsdifferenzenverfahren fällt linear mit kleiner werdender Schrittweite h. Für asymptotische Funktionsverläufe können je nach gewünschter Genauigkeit etwa 100 bis 1000 Stützstellen nötig sein. Schwingende Systeme mit Differenzialgleichungen 2. Ordnung brauchen oft deutlich kleinere Schrittweiten, besonders wenn mehrere Perioden einer gedämpften Schwingung dargestellt werden sollen. Mit kleiner werdendem h gleichen sich Vorwärts- und Rückwärtsverfahren bei gleichem Eingangssignal an. Der Gesamtfehler besteht aus dem Fehler des numerischen Verfahrens und dem Rundungsfehler; ungenaue Zahlenwerte können Rundungsfehler über viele Folgeglieder aufsummieren. Laut einer zitierten Veröffentlichung von Thomas Westermann ergab eine Fehleranalyse bei h=5·10^{-5} s für Euler-Vorwärts gegenüber Runge-Kutta etwa den 12-fachen absoluten Fehler.

Nichtlineare Übertragungssysteme sind mit konventionellen Methoden im kontinuierlichen Zeitbereich oft nicht geschlossen lösbar. Mit Differenzengleichungen können sie näherungsweise berechnet werden. Statische Nichtlinearitäten lassen sich zum Beispiel durch IF-THEN-ELSE-Anweisungen, Tabellen, Begrenzungen oder gespeicherte Hysteresefunktionen beschreiben. Nach dem Hammersteinmodell wird ein statisches nichtlineares Modell in Reihe mit einem linearen dynamischen System betrachtet.

Ein weiteres Sonderthema ist die Systemtotzeit. Sie kann durch Annäherungsmodelle mit Verzögerungen höherer Ordnung oder Allpassgliedern behandelt werden. Numerisch ist auch die Speicherung von Folgegliedern möglich, sodass auf Werte zugegriffen wird, die um die Totzeit T_t zurückliegen. In Tabellenkalkulationen kann dafür etwa eine INDEX-Funktion verwendet werden. Für Simulationen dynamischer Übertragungssysteme sind Testsignale wichtig, zum Beispiel Sprungfunktion, Impulsfunktion oder Anstiegsfunktion. Diese nichtperiodischen deterministischen Testsignale beginnen bei t=0 und haben für t<0 den Wert u(t) ≡ u_k = 0.