vektorwertiges DGL-System 2. Ordnung mit Matlab lösen

WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

saint

Forum-Newbie


	Beiträge: 6

	Anmeldedatum: 18.09.09

	Wohnort: ---

	Version: ---

Verfasst am: 23.05.2011, 11:47 Titel: vektorwertiges DGL-System 2. Ordnung mit Matlab lösen

Hallo,

ich experimentiere nun schon seit einigen Tagen in Matlab, um folgendes vektorwertige Differentialgleichungssystem (nachfolgend kurz DGLS) mit

zu lösen:

Ich weiß, dass man das DGLS in ein DGLS erster Ordnung umwandeln muss. Dies erscheint mir hier auch sonderbar einfach:

Wie bekomme ich nun allerdings dieses vektorwertige DGLS erster Ordnung mit Matlab gelöst? Bislang habe ich folgenden Code geschrieben, der aber definitiv nicht in Ordnung ist (Matlab rechnet > 5 Minuten):

Code:

function xr = diffEq(t,x)

mu0 = 1.32712440018000e+20;
mu1 = 2.203200000000000e+013;
mu2 = 3.248590000000000e+014;

xr = zeros(size(x));

hx0 = [x(1);x(2);x(3)];
hx1 = [x(7);x(8);x(9)];
hx2 = [x(13);x(14);x(15)];

ddr = f(mu1,hx0,hx1) + f(mu2,hx0,hx2);
xr(1) = ddr(1);
xr(2) = ddr(2);
xr(3) = ddr(3);
ddr = f(mu0,hx1,hx0) + f(mu2,hx1,hx2);
xr(4) = ddr(1);
xr(5) = ddr(2);
xr(6) = ddr(3);
ddr = f(mu0,hx2,hx0) +f(mu1,hx2,hx1);
xr(7) = ddr(1);
xr(8) = ddr(2);
xr(9) = ddr(3);
end

function out = f(mu,pm,pM)
% returns a column vector out = [x;y;z]
out = -((mu/(norm(pm-pM)))*(pm-pM));
end

Funktion ohne Link?

Der von Matlab übergebene Parameter x hat die Struktur (Spalten) $r_{0,x}$ , $r_{0,y}$ , $r_{0,z}$ , $\dot{r}_{0,x}$ , $\dot{r}_{0,y}$ , $\dot{r}_{0,z}$ , $r_{1,x}$ , $r_{1,y}$ , $r_{1,z}$ , $\dot{r}_{1,x}$ , $\dot{r}_{1,y}$ , $\dot{r}_{1,z}$ , $r_{2,x}$ , $r_{2,y}$ , $r_{2,z}$ , $\dot{r}_{2,x}$ , $\dot{r}_{2,y}$ , $\dot{r}_{2,z}$ . Aus diesem Grunde habe ich immer drei Spalten x,y,z zu einem Vektor zusammengefasst (hx0 für r0, hx1 für r1, hx2 für r2).

Ich scheine allerdings Probleme dabei zu haben, den Output-Vektor für ode45, xr, zusammenzusetzen und vermute dort meinen Fehler.

Über jede Hilfe würde ich mich freuen.

Beste Grüße,
saint

P.S.: In Simulink konnte ich dieses DGLS bereits lösen. Ich wollte eigentlich mit der Lösung über ein Matlab-Skript erreichen, dass man nicht so viele Verdrahtungen in Simulink durchführen muss.

Harald

Forum-Meister


	Beiträge: 24.492

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 23.05.2011, 16:39 Titel:

Hallo,

stimmt denn die mit MATLAB alleine berechnete Lösung am Ende überhaupt? Ich vermute nicht.

Problem: du belegst die Komponenten 10-18 von xr nicht.

Meines Erachtens müsste es so aussehen:

Code:

ddr = f(mu1,hx0,hx1) + f(mu2,hx0,hx2);
xr(1:3) = x(4:6)';
xr(4) = ddr(1);
xr(5) = ddr(2);
xr(6) = ddr(3);

Funktion ohne Link?

etc.

Grüße,
Harald

saint

Themenstarter

Forum-Newbie


	Beiträge: 6

	Anmeldedatum: 18.09.09

	Wohnort: ---

	Version: ---

Verfasst am: 23.05.2011, 17:45 Titel:

Hallo Harald,

vielen Dank für deine rasche Antwort. Ich kann leider nicht sagen, ob die berechnete Lösung stimmt oder nicht, denn ich habe die Berechnung nach 10 Minuten abgebrochen und habe somit keine Ergebnisse (in Simulink dauerte die Berechnung nicht einmal eine Sekunde).

Ich habe den Code nun wie folgt geändert (entsprechend deines Hinweises):

Code:

function xr = diffEq(t,x)

mu0 = 1.32712440018000e+20;
mu1 = 2.203200000000000e+013;
mu2 = 3.248590000000000e+014;

xr = zeros(size(x));

hx0 = [x(1);x(2);x(3)];
hx1 = [x(7);x(8);x(9)];
hx2 = [x(13);x(14);x(15)];

ddr = f(mu1,hx0,hx1) + f(mu2,hx0,hx2);
xr(1) = x(4);
xr(2) = x(5);
xr(3) = x(6);
xr(4) = ddr(1);
xr(5) = ddr(2);
xr(6) = ddr(3);
ddr = f(mu0,hx1,hx0) + f(mu2,hx1,hx2);
xr(7) = x(10);
xr(8) = x(11);
xr(9) = x(12);
xr(10) = ddr(1);
xr(11) = ddr(2);
xr(12) = ddr(3);
ddr = f(mu0,hx2,hx0) +f(mu1,hx2,hx1);
xr(13) = x(16);
xr(14) = x(17);
xr(15) = x(18);
xr(16) = ddr(1);
xr(17) = ddr(2);
xr(18) = ddr(3);
end

function out = f(mu,pm,pM)
% returns a column vector out = [x;y;z]
out = -((mu/(norm(pm-pM)))*(pm-pM));
end

Funktion ohne Link?

Aufgerufen wird ode45 über

Code:

[t,x] = ode45(@diffEq, [0:86400:86400*2], [r0' dr0' r1' dr1' r2' dr2']);

Funktion ohne Link?

.

Leider läuft auch dieser Versuch nun schon seit über fünf Minuten ohne Ergebnis. Hast du noch eine Idee, worin das Problem liegen könnte?

Beste Grüße,
saint

Harald

Forum-Meister


	Beiträge: 24.492

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 23.05.2011, 18:17 Titel:

Hallo,

gib doch bitte die Anfangsbedingungen an, damit man es hier selbst ausprobieren kann.

Wo ist der Exponent 3 aus deinen Gleichungen in deiner Funktion abgeblieben?
Hast du in Simulink etwas an den Einstellungen geändert?
Hast du mal eine kurze Simulation gemacht und die Ergebnisse verglichen?

Im übrigen könnte man den Code sicher effizienter gestalten:

Code:

x(1:3) = ...
x(4:6) = ...

Funktion ohne Link?

Grüße,
Harald

saint

Themenstarter

Forum-Newbie


	Beiträge: 6

	Anmeldedatum: 18.09.09

	Wohnort: ---

	Version: ---

Verfasst am: 23.05.2011, 18:42 Titel:

Hallo Harald,

abermals vielen Dank für deine Hilfe. In der Tat hatte ich den Exponenten vollkommen übersehen - und das war auch die Fehlerursache! Nachdem die Funktion f(..) nachkorrigiert ist, funktioniert nun alles wie gewünscht und ich bekomme korrekte Ergebnisse! Very Happy

Ich werde den Code noch kompakter schreiben; zum ersten Ausprobieren wollte ich allerdings alle Dinge ausführlich schreiben, so dass ich im Nachhinein noch alles kommentieren kann. Bislang war ich noch nicht mit der Lösung von DGLs in Matlab in Berührung gekommen.

Für deine Hilfe danke ich dir noch einmal recht herzlich und wünsch dir alles Gute!

saint


Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:	Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen

Impressum | Nutzungsbedingungen | Datenschutz | FAQ |

RSS

Hosted by:

Copyright © 2007 - 2024 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.