Initial value problems for system of differential-algebraic equations in Maple

Anwendungen von DAEs treten natürlich in vielen Bereichen auf, z.B. verschiedene dynamische Prozesse, mechanische Systeme, Simulation von elektrischen Schaltkreisen und chemischen Reaktionen, die Invarianten unterliegen usw., und diese werden oft durch DAEs ausgedrückt, die aus algebraischen Gleichungen und Differentialoperationen bestehen. Mehrere Methoden und Algorithmen sind von vielen Forschern und Ingenieuren eingeführt worden, um die IVPs für Systeme von DAEs zu lösen. Die meisten von ihnen versuchen, eine ungefähre Lösung des gegebenen Systems zu finden. Wir erinnern jedoch an einen symbolischen Algorithmus, um die exakte Lösung eines gegebenen Systems von DAEs zu berechnen (siehe für weitere Details des Algorithmus). In diesem Beitrag besprechen wir das Maple-Paket des symbolischen Algorithmus, der die exakte Lösung berechnet.

Es gibt mehrere implementierte Methoden, die in verschiedenen mathematischen Softwares wie Matlab, Mathematica, SCIlab usw. verfügbar sind. Alle diese Implementierungen werden angewendet, um die allgemeine Lösung eines gegebenen Systems DAEs mit freien Parametern zu finden. Dann können wir die Werte der Parameter durch die Substitution der Anfangsbedingungen finden. Die in Mathematica implementierte Methode basiert zum Beispiel auf der Zerlegung der Koeffizientenmatrizen A und B in einen nonsingulären und einen nilpotenten Teil. Dann wird die verallgemeinerte Inverse für A und B berechnet, und das Problem wird auf die Lösung eines Systems von ODEs reduziert. Es können also vorhandene Löser für ODEs verwendet werden. In Matlab wird die Gleichung auch in ein System von ODEs umgewandelt, indem der Differentialindex reduziert wird, und dann finden wir die allgemeine Lösung mit freien Parametern. In dem vorgeschlagenen Algorithmus berechnen wir jedoch die exakte Lösung direkt ohne freie Parameter. Das implementierte Maple-Paket basiert auf der Umwandlung des gegebenen Systems in eine kanonische Form mit Hilfe des Shuffle-Algorithmus, der ein weiteres einfaches äquivalentes System erzeugt, und das kanonische System kann leicht gelöst werden. Die Maple-Implementierung beinhaltet die Berechnung des kanonischen Systems und der exakten Lösung eines gegebenen IVP. Der Vergleich des implementierten Maple-Pakets mit bestehenden Methoden, die in anderer mathematischer Software wie Matlab und Mathematica implementiert sind, wird auch im Abschnitt Ergebnisse diskutiert. In dieser Arbeit haben wir uns auf die Maple-Implementierung eines IVP mit homogenen Anfangsbedingungen konzentriert, wir diskutieren jedoch auch einen Algorithmus zur Überprüfung der Konsistenz der nicht-homogenen Anfangsbedingungen.

Symbolischer Algorithmus von IVPs für Systeme von lDAEs

In dieser Arbeit haben wir uns auf ein System von DAEs der allgemeinen Form

$$\begin{aligned} A{\texttt{D}}y(x)+By(x)=f(x). \end{aligned}$

(1)

Das System (1) ist ein rein algebraisches System, wenn $A = 0$, und es gibt viele Methoden und Algorithmen, um alle möglichen Lösungen zu berechnen, siehe z.B. . Das System (1) wird zu einem System von gewöhnlichen Differentialgleichungen, wenn $\det (A) \ne 0$ (wir nennen A eine reguläre Matrix), und die Lösung eines Systems gewöhnlicher Differentialgleichungen erster Ordnung wird in diskutiert. Daher haben wir uns auf ein System der Form (1) konzentriert, bei dem die Koeffizientenmatrix A eine singuläre Matrix ungleich Null ist. Nehmen wir der Einfachheit halber an, dass $\mathcal {F} = C^\infty $ und $ \subset \mathbb {R}$. Nun kann die Operatorform eines IVP für DAEs wie folgt dargestellt werden

$\begin{aligned}&Ly = f, \\&{\texttt{E}}y = 0, \end{aligned}$

wobei $L = A{\texttt{D}}+ B \in \mathcal {F}^{n \times n}$ der Matrixdifferentialoperator ist, ${\texttt{D}}= \frac{d}{dx}$, $y \in \mathcal {F}^n$ ist der zu bestimmende unbekannte Vektor, $f \in \mathcal {F}^n$ ist die Vektorzwangsfunktion und $\texttt{E}$ ist der Auswerteoperator. Wir wollen den Matrix-Green-Operator $G \in \mathcal {F}^{n \times n}$ so finden, dass $Gf = y$ und ${\texttt{E}}G = 0$.

Das folgende Lemma 1 ist einer der wesentlichen Schritte für den vorgeschlagenen Algorithmus. Das Lemma gibt die Variation der Parameterformel eines IVP für skalare lineare Differentialgleichungen höherer Ordnung über Integro-Differential-Algebren an.

Lemma 1

(Thota und Kumar ) Sei $(\mathcal {F}, {\texttt{D}}, {\texttt{A}})$ eine gewöhnliche Integro-Differential-Algebra. Nehmen wir an, $T = {\texttt{D}}^m + a_{m-1}{\texttt{D}}^{m-1} + \cdots + a_0 \in \mathcal {F}$ sei ein monischer skalarer Differentialoperator der Ordnung m und $v_1, \ldots , v_m$ sei ein Fundamentalsystem für T. Dann ist der rechtsinverse Operator von T gegeben durch

$\begin{aligned} T^\divideontimes = \sum \limits _{i=1}^{m} v_i {\texttt{A}}w^{-1}w_i \in \mathcal {F}, \end{aligned}$

(2)

wobei w die Determinante der Wronskischen Matrix W für $v_1, \, v_{m}$ und $w_i$ die Determinante der Matrix $W_i$, die sich aus W durch Ersetzen der i-ten Spalte durch den m-ten Einheitsvektor ergibt.

Um die Green’sche Funktion und die exakte Lösung eines gegebenen DAE-Systems zu erhalten, finden wir zunächst eine kanonische Form des gegebenen DAE-Systems unter Verwendung des Shuffle-Algorithmus, der das gegebene System in ein anderes äquivalentes und einfacheres System transformiert, das leicht gelöst werden kann. Der Green’sche Operator und die Green’sche Funktion des gegebenen Systems von DAEs mit Anfangsbedingungen werden mit Hilfe der kanonischen Form im folgenden Theorem berechnet.

Theorem 2

(Thota und Kumar ) Sei ($\mathcal {F}, {\texttt{D}}, {\texttt{A}})) eine gewöhnliche Integro-Differential-Algebra. Sei \(\tilde{L} = \tilde{A}{\texttt{D}}+ \tilde{B} \in \mathcal {F}^{n \times n}$ die kanonische Form von $L = A{\texttt{D}}+B$ mit Anfangsbedingungen; und $\{v_1, \ldots , v_{n}\}$ sei ein fundamentales System für $T = \det (\tilde{L})$. Dann ist das reguläre IVP für das System der DAEs

$\begin{aligned}&Ly = f, \\&{\texttt{E}}y = 0, \end{aligned}$

hat die eindeutige Lösung

$$\begin{aligned} y = \begin{pmatrix} \sum _{i=1}^{n} (-1)^{i+1} \mathcal {L}_i^1 T^\divideontimes \tilde{f_i} \\ (-1)^{i+1} \mathcal {L}_i^1 (-1)^{i+n} \mathcal {L}_i^n T^\divideontimes \tilde{f_i} \end{pmatrix}, \end{aligned}$

(3)

wobei $\mathcal {L}_i^j$ die Determinante von $\tilde{L}$ nach Entfernen der i-ten Zeile und j-ten Spalte ist; $T^\divideontimes$ ist die rechte Inverse von T; und $\tilde{f} = (\tilde{f_1}, \ldots , \tilde{f_n})^t$. Der Greensche Operator ist

$\begin{aligned} G = \begin{pmatrix} (-1)^{1+1} \mathcal {L}_1^1 T^\divideontimes &{} \cdots &{} (-1)^{n+1} \mathcal {L}_n^1 T^\divideontimes \\ \vdots &{} \ddots &{} \vdots \\(-1)^{1+n} \mathcal {L}_1^n T^\divideontimes &{} \cdots &{} (-1)^{n+n} \mathcal {L}_n^n T^\divideontimes \end{pmatrix} \end{aligned}$

(4)

so dass $G\tilde{f} = y$ und ${\texttt{E}}~G = 0$.

Nicht-homogene Anfangsbedingungen

Im Allgemeinen gibt es keine Freiheit, nicht-homogene Anfangsbedingungen für die vorgeschlagene Methode zu wählen. Daher haben die Autoren einen Algorithmus vorgestellt, um die Konsistenz einer gegebenen inhomogenen Anfangsbedingung zu überprüfen. In diesem Abschnitt rufen wir den Algorithmus in Proposition 3 auf, um die Konsistenz der nicht-homogenen Anfangsbedingungen zu prüfen.

Proposition 3

(Thota und Kumar ) Sei ($\mathcal {F}, {\texttt{D}}, {\texttt{A}}$) eine gewöhnliche Integro-Differential-Algebra. Nehmen wir an, $\tilde{T} = \tilde{A}{\texttt{D}}+ \tilde{B} \in \mathcal {F}^{n \times n}$ ist eine kanonische Form von $T = A{\texttt{D}}+B$ mit nicht-homogenen Anfangsbedingungen. Die inhomogene Anfangsbedingung ${\texttt{E}}u = \alpha$ ist konsistent, wenn

$$begin{aligned} UU_a^{-1}\alpha \in \text {Ker}(T), \end{aligned}$

(5)

Wobei U die Fundamentalmatrix von $\tilde{T}$ ist und $U_a$ der Wert von U am Anfangspunkt a ist.

Das folgende Beispiel zeigt die Berechnung der exakten Lösung unter Verwendung des in Satz 2 vorgestellten Algorithmus und prüft außerdem die Konsistenz der inhomogenen Anfangsbedingungen unter Verwendung des in Satz 3 vorgestellten Algorithmus.

Beispiel 4

Betrachten Sie die folgenden differential-algebraischen Gleichungen.

$$begin{aligned}y_1’+y_2’+y_1+y_2&=x \\y_1-y_2 &= \sin x \end{aligned}$

(6)

mit Anfangsbedingung $y_1(0)=y_2(0)=0$.

Der Matrixdifferentialoperator und die kanonische Form von (6) sind

$$begin{aligned}T = \begin{pmatrix} 1+{\texttt{D}}~&{}~ 1+{\texttt{D}}\\ 1 &{} -1 \end{pmatrix},~~~ \tilde{T} = \begin{pmatrix} {\texttt{D}}&{} 2+{\texttt{D}}\ {\texttt{D}}&{} -{\texttt{D}}\end{pmatrix} \text {und}~ \tilde{f} = \begin{pmatrix} x-\sin x \\\ \cos x \end{pmatrix}. \end{aligned}$

Nach dem Algorithmus in Theorem 2 haben wir die exakte Lösung

$$begin{aligned} y = \begin{pmatrix} \frac{1}{2}e^{-x} + \frac{1}{2} \sin x + \frac{1}{2} x – \frac{1}{2}\ \frac{1}{2}e^{-x} – \frac{1}{2} \sin x + \frac{1}{2} x – \frac{1}{2} \end{pmatrix}. \end{aligned}$

(7)

Man kann leicht überprüfen, dass $Ty = f$ und ${\texttt{E}}y = 0$.

Betrachten Sie die inhomogenen Anfangsbedingungen $y_1(0) = \alpha , y_2(0)=\beta$ mit gegebenem System (6). Aus Proposition 3 sind die Anfangsbedingungen konsistent, wenn $UU_0^{-1} \alpha \in \text {Ker}(T)$, wobei

$$\begin{aligned} U=\begin{pmatrix} 1 &{} e^{-x} \\ 0 &{} e^{-x} \end{pmatrix}, U_0= \begin{pmatrix}1 &&{} 1 \end{pmatrix},~~\text {and}~~~ \alpha = \begin{pmatrix}\alpha \\\ \beta \end{pmatrix}. \end{aligned}$

Nun

$$\begin{aligned} UU_0^{-1} \alpha = \begin{pmatrix} 1 &{} e^{-x} \\ 0 &{} e^{-x} \end{pmatrix} \begin{pmatrix}1 &&{} 1 \end{pmatrix} \begin{pmatrix}\alpha \\ \beta \end{pmatrix} = \begin{pmatrix} e^{-x}\beta + \alpha – \beta \\\ e^{-x}\beta \end{pmatrix} \end{aligned}$

und

$$\begin{aligned} T(UU_0^{-1} \alpha ) = \begin{pmatrix} 1+{\texttt{D}}~&{}~ 1+{\texttt{D}}\\ 1 &{} -1 \end{pmatrix} \begin{pmatrix} e^{-x}\beta + \alpha – \beta \\\ e^{-x}\beta \end{pmatrix} = \begin{pmatrix} \alpha – \beta \\alpha – \beta \end{pmatrix}, \end{aligned}$

Das ergibt $\alpha – \beta = 0$ für $UU_0^{-1} \alpha \in \text {Ker}(T)$, und daher sind die konsistenten Anfangsbedingungen $y_1(0) = \alpha ,y_2(0)=\beta$, so dass $\alpha -\beta =0$. Die Lösung $y_c$ des IVP $Ty=0, {\texttt{E}}y = (\alpha ,\beta )^T$ berechnet sich als (siehe für weitere Details),

$$\begin{aligned} y_c = UU_0^{-1} \alpha = \begin{pmatrix} e^{-x}\beta + \alpha – \beta \\\ e^{-x}\beta \end{pmatrix}, \end{aligned}$

und aus (7)

$$\begin{aligned}y_p=\begin{pmatrix} \frac{1}{2}e^{-x} + \frac{1}{2} \sin x + \frac{1}{2} x – \frac{1}{2}\ \frac{1}{2}e^{-x} – \frac{1}{2} \sin x + \frac{1}{2} x – \frac{1}{2} \end{pmatrix}. \end{aligned}$

Die exakte Lösung des regulären Systems

$$\begin{aligned} \{\begin{array}{*{20}c} {1 + {\texttt{D}}} & {1 + {\texttt{D}} \\\ 1 & { – 1} \\end{array} } \\end{array} } \right) & = \left( {\begin{array}{*{20}c} {y_{1} } \\\\ {y_{2} } \\\end{array} } \right), \\ E\left( {\begin{array}{*{20}c} {y_{1} } \\\ {y_{2} } \\\ end{array} } \right) & = \left( {\begin{array}{*{20}c} \alpha \\\ \beta \\\ \end{array} } \right),\,\,\alpha – \beta = 0, \\\ \end{aligned}$

ist $y=y_c+y_p$, d.e.,

$$begin{aligned} y = \begin{pmatrix} \frac{1}{2}e^{-x} + \frac{1}{2} \sin x + \frac{1}{2} x – \frac{1}{2}+e^{-x}\beta + \alpha – \beta \\\ \frac{1}{2}e^{-x} – \frac{1}{2} \sin x + \frac{1}{2} x – \frac{1}{2}+e^{-x}\beta \end{pmatrix}. \end{aligned}$

Heading