DAE aus MNA lösen (nur am Rande Programmierung, aber hier hängen die Mathematiker rum

[zgep]

Hi zusammen!

Gleich mal als Entschuldigung: Zwar ist das Problem eher ein mathematisches, gelöst wird es aber mittels Matlab-, Fortran- und C++-Code ;)

Ich arbeite seit Mo an einem neuen Projekt mit (icestars.eu, WP1 - Time Domain). Als eher FPGA-Designer und embedded Programmierer sind da für mich ziemlich viele neue Konzepte dabei.

Derzeit steige ich bei den Methoden zur Lösung der DAEs nicht ganz durch. Ich würde das ganze aber ganz gern verstehen, anstatt nur die fertigen Lösungen zu verwenden.


Grundsätzlich sieht das Problem wie folgt aus:

f(x(t)) + q'(x(t)) = 0, wobei f, x, q Element aus R^n

f, q, x(t_n), t_n, delta_t sind bekannt, gesucht ist x(t_n + delta_t) = x(t_n) + delta_x


Die fertige Lösung sieht so aus:

delta_x = -(Jf + Jq/delta_t)\(f + (q - q0)/delta_t)

Jf und Jq sind die Jacobi-Matrizen zu f und q. Natürlich ist es mit einer Iteration nicht getan, aber das ist nicht der Teil den ich nicht verstehe.

falls \ unbekannt ist: in Matlab bedeutet A\b das selbe wie inv(A)*b.


Mit etwas umformen und Hilfe von Herrn Euler ergibt sich

-delta_x*(Jf + Jq/delta_t) = f + q'

ausmultipliziert und f auf beiden Seiten weggekürzt ist das dann

-delta_x*Jq/delta_t = q'


ich weis, dass -Jq*delta_x' = q' und vermute, dass man daraus irgendwie darauf schließen kann, dass das obige korrekt ist. Aber wie das funktioniert verstehe ich nicht.

Kann mir dabei jemand weiterhelfen?

Danke!
zgep

[Pinoccio]

falls \ unbekannt ist: in Matlab bedeutet A\b das selbe wie inv(A)*b.sOT: VOrsicht, das kann schiefgehen (http://www.tm-mathe.de/Themen/html/matlab__probleme_mit_singulare.html).

Ansonsten, DNA hab ich nur mal angerissen (in TP), aber auch nicht numerisch.

mfg

[Gast]

Das die Matrix ausreichend konditioniert ist wird vom C++ Netzwerkanalysator sichergestellt.
Aber das der Hinweis überhaut kommt ist schon mal ein gutes Zeichen, dass ich hier die richtigen Leute befrage!

mfg,
zgep

[Pinoccio]

Das die Matrix ausreichend konditioniert ist wird vom C++ Netzwerkanalysator sichergestelltGut. I. d. R. ist das eh kein Problem, singuläre Matritzen sind zwar u. U. ein numerischer Alptraum, wenn man sie vergisst, oftmals entspringen sie aber einfach nur einem schlechten Design.

(1) -delta_x*Jq/delta_t = q'
(2) -Jq*delta_x' = q'
und vermute, dass man daraus irgendwie darauf schließen kann, dass das obige korrekt ist. Aber wie das funktioniert verstehe ich nicht.

Kann mir dabei jemand weiterhelfen?Falls deine Frage lautet, wie man aus (1) auf (2) kommt, dann ist die Antwort womöglich
delta_x'=delta_x/delta_t
der 1/8-Physiker in mir schreit dann aber, daß man Zeitableitungen punktet.

mfg

[zgep]

jepp, das ist die Frage. Sorry wegen der Ticks, aber ich habe keine Ahnung wie ich Punkte drüber machen könnte und jedesmal d(f)/dt zu schreiben war mir zu aufwendig ;)
leider ist delta_x' eben nicht gleich delta_x/delta_t (zumindest sehe ich nicht wie das zutreffen kann). Es muss irgendwas mit der Linearisierung durch die Jacobi-Matrix zu tun haben. Genau da steig ich leider nicht durch. delta_x/delta_t entspricht ja eigentlich G' und delta_x' G'' ?

Ich kann mir zwar die allgemeine Formel wie sie auf der wikipedia steht

F(p + delta) = F(p) + JF*delta

, daraus ergibt sich zB das delta_x = -Jf\f bei der OP-Berechnung im selben Projekt wo der q-Vektor 0 ist ) erklären, aber warum

-Jq*delta_x/delta_t = d(q)/dt

sein soll versteh ich nicht.

Dennoch natürlich danke für die Antworten!
zgep

[Gast]

Nachtrag:

etwas anders formuliert steht da ja

-Jq*( x(t + delta_t) - x(t) )/delta_t = q' = -Jq*d( x(t + delta_t) - x(t) )/dt

( x(t + delta_t) - x(t) )/delta_t = d( x(t + delta_t) - x(t) )/dt

und ich sehe nicht wie das gleich sein könnte.

mfg,
zgep

[Pinoccio]

jepp, das ist die Frage. Sorry wegen der Ticks, aber ich habe keine Ahnung wie ich Punkte drüber machen könnte und jedesmal d(f)/dt zu schreiben war mir zu aufwendig ;)Achso. :biggrin:
Falls du Latex beherrschst, kannst du hier (http://mathbin.net/) deine Formeln (auch mehrzeilig) eintippen und das fertige Bild hier mit [ img]...[/img] einbinden.
Beispiel:
http://mathbin.net/equations/12355_0.png (<-link) (http://mathbin.net/12355)
leider ist delta_x' eben nicht gleich delta_x/delta_t (zumindest sehe ich nicht wie das zutreffen kann). [...]Das vertehe ich nicht. Was ist denn delta_x' dann? Die Ableitung von x, von delta_x und wonach oder garkeine Ableitung?

Dennoch natürlich danke für die Antworten!Kein Problem.

/edit:
Nachtrag:

etwas anders formuliert steht da ja

-Jq*( x(t + delta_t) - x(t) )/delta_t = q' = -Jq*d( x(t + delta_t) - x(t) )/dt

( x(t + delta_t) - x(t) )/delta_t = d( x(t + delta_t) - x(t) )/dt

und ich sehe nicht wie das gleich sein könnte.

mfg,
zgepLinks und rechts der unteren Gleichung?, also quasi http://mathbin.net/equations/43979_0.png? (http://mathbin.net/43973)
/edit: Bild gefixt


mfg

[zgep]

http://mathbin.net/equations/12355_0.png (http://mathbin.net/12355) (<-link)*Bookmark*Was ist denn delta_x' dann? Die Ableitung von x, von delta_x und wonach oder garkeine Ableitung?

delta_x ist x(t + delta_t) - x(t). Dadurch ist delta_x/delta_t näherungsweise die Ableitung von x(t) nach t.

delta_x' ist die Ableitung von delta_x nach t.


es steht da quasi (wenn man die J-Matrix einfach wegkürzt)

(ich bin zu blöd das Bild zu finden) (http://mathbin.net/43975)

mfg,
zgep

[zgep]

/edit:
Links und rechts der unteren Gleichung?, also quasi http://mathbin.net/equations/43973_0.png (http://mathbin.net/43973)?

genau! nur ne Klammer fehlt ;)

mfg,
zgep

[Pinoccio]

delta_x ist x(t + delta_t) - x(t). Dadurch ist delta_x/delta_t näherungsweise die Ableitung von x(t) nach t.

delta_x' ist die Ableitung von delta_x nach t.


es steht da quasi (wenn man die J-Matrix einfach wegkürzt)Hm. Vielleicht fällt mir morgen dazu noch was ein, heute seh ich nichts mehr.
Vielleicht kommt ja auch noch pest vorbei oder ein anderer.
(ich bin zu blöd das Bild zu finden) (http://mathbin.net/43975)Entweder über Rechtsklick und Bildeigenschaften oder
[ url=http://mathbin.net/43975][ img]http://mathbin.net/equations/43975_0.png[/img][/url] /edit, stimmt nur für das erste Bild, bei replys nicht.

/edit:
genau! nur ne Klammer fehlt ;)gefixt (http://mathbin.net/43973)

mfg