Lösung von DGL mit ode45

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

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

yiiit

Forum-Newbie


	Beiträge: 8

	Anmeldedatum: 21.10.11

	Wohnort: Osnabrück

	Version: ---

Verfasst am: 01.12.2011, 11:46 Titel: Lösung von DGL mit ode45

Hallo!

Ich möchte die folgenden DGLs mit der Matlab-Funktion ode45 lösen, habe aber keine Ahnung, wie man die zugehörigen odefunktionen aufstellt. Die DGLs lauten:

$\dot{v_{x}} = g \cdot \frac{r_{E}^{2}}{(x^{2}+y^{2})^{\frac{3}{2}}} \cdot x$

$\dot{x} = v_{x}$

$\dot{v_{y}} = g \cdot \frac{r_{E}^{2}}{(x^{2}+y^{2})^{\frac{3}{2}}} \cdot y$

$\dot{y} = v_{y}$

offensichtlich muss man zwei odefunktionen schreiben, eine für die ersten beiden gleichungen und eine für die letzten beiden gleichungen - aber ich hab keine ahnung wie das funktioniert.

danke!

Harald

Forum-Meister


	Beiträge: 24.502

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 01.12.2011, 12:40 Titel:

Hallo,

nein, du kannst mit einer Funktion arbeiten.

In der Dokumentation ist ein Beispiel, das deinem von der Struktur her recht ähnlich ist:
http://www.mathworks.com/help/relea.....b/techdoc/ref/ode23t.html
oder

Code:

doc ode23

Funktion ohne Link?

(die Infos gelten auch für ode45)

Grüße,
Harald

yiiit

Themenstarter

Forum-Newbie


	Beiträge: 8

	Anmeldedatum: 21.10.11

	Wohnort: Osnabrück

	Version: ---

Verfasst am: 01.12.2011, 12:56 Titel:

offensichtlich handelt es sich bei meinen Gleichungen ja um einen Vektor, den ich zerlegt habe aus einer DGL 2. Ordnung.. Danke für deine schnelle Hilfe, Harald, aber ich finde bei den Beispielen nichts Ähnliches, sprich mit drei Variablen.. ich brauche aber bei der ode-Funktion doch die Variablen t, y, x - ich finde aber immer nur y und t.

Jan S

Moderator


	Beiträge: 11.057

	Anmeldedatum: 08.07.10

	Wohnort: Heidelberg

	Version: 2009a, 2016b

Verfasst am: 01.12.2011, 13:23 Titel:

Hallo yiiit,

Dann sind die beiden Komponenten x und y zu einem Vektor zusammengefasst, der in der Dokumentation y heißt. Man kann dann leicht ein DLG-System 1. Ordnung in den Variablen y(1) (also Dein "y") und y(2) (Dein "x", oder umgekehrt) und der Zeit t erstellen.

Gruß, Jan

yiiit

Themenstarter

Forum-Newbie


	Beiträge: 8

	Anmeldedatum: 21.10.11

	Wohnort: Osnabrück

	Version: ---

Verfasst am: 01.12.2011, 13:43 Titel:

also:

Code:

function dot = DGL(t,y)

Re = 6371;
g = 9.81/1000;

dot = [0;0;0;0];

dot(1) = y(3);
dot(2) = y(4);
dot(3) = g.*Re.^2./((y(1).^2+y(2).^2).^(3/2)).*y(1);
dot(4) = g.*Re.^2./((y(1).^2+y(2).^2).^(3/2)).*y(1);

Funktion ohne Link?

und

Code:

function [] = main

[T,Y] = ode45(@DGL,[0 20000], [10000 0 -8*cos(0) 8*sin(0)])

plot(T,Y)

Funktion ohne Link?

Aber wie plotte ich das denn jetzt so, dass eine Bahnkurve existiert? ich hab ja jetzt nur wild durch den Raum fliegende Geraden und Linien..

Jan S

Moderator


	Beiträge: 11.057

	Anmeldedatum: 08.07.10

	Wohnort: Heidelberg

	Version: 2009a, 2016b

Verfasst am: 01.12.2011, 13:51 Titel:

Hallo yiiit,

Was genau sind "wild durch den Raum fliegende Geraden und Linien"?

Möglicherweise benötigst Du noch Zwischenwerte, die Du leicht im ODE45-Input TSPAN angeben kannst.

Gruß, Jan

Harald

Forum-Meister


	Beiträge: 24.502

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 01.12.2011, 14:03 Titel:

Hallo,

zusätzlich zu Jan's Vorschlag will man bei einer Bahnkurve normal nicht die Komponenten in Abhängigkeit der Zeit plotten, sondern x- und y-Werte gegeneinander:

Code:

plot(Y(:,1), Y(:,2)) % ... oder animiert
comet(Y(:,1), Y(:,2))

Funktion ohne Link?

Grüße,
Harald

yiiit

Themenstarter

Forum-Newbie


	Beiträge: 8

	Anmeldedatum: 21.10.11

	Wohnort: Osnabrück

	Version: ---

Verfasst am: 01.12.2011, 14:14 Titel:

Also der Hintergrund der Aufgabe ist ja, dass wir für einen bestimmten Winkel alpha die Trajektorie der Bahnkurve einer zur Erde zurückkommenden Raumkapsel bestimmen sollen..

Ich hab eigentlich grad Probleme damit, mein Ergebnis zu interpretieren.. Eigentlich will ich ja einen Vektor r, der mir für jeden Zeitpunkt die x- und die y-Koordinate ausspuckt mitsamt der zugehörigen Geschwindigkeit..

Steht jetzt in meinem y-Vektor an der Stelle y(1) der x-Wert und an der Stelle y(2) der y-Wert, an der Stelle y(3) meine Geschwindigkeit in x-Richtung und an der Stelle y(4) meine Geschwindigkeit in y-Richtung?

Aber wenn ich in den Anfangsbedingungen einen Winkel von 0° eintrage, dann müsste doch die Kapsel gerade auf die Erde fliegen und wenn ich mir den Plot von x gegen y dann anschaue, sehe ich, dass die Kapsel kurz vorher so stark abgelenkt wird und umdreht?

Ich brauche also Hilfe beim Interpretieren des Ergebnisses und es wäre nett, wenn du mir nochmal deinen Tipp mit den Zwischenwerten genauer erklären könntest.

Aber meine DGL.m und meine main.m sind sonst soweit den eingangs beschriebenen DGLs entsprechend korrekt?

Jan S

Moderator


	Beiträge: 11.057

	Anmeldedatum: 08.07.10

	Wohnort: Heidelberg

	Version: 2009a, 2016b

Verfasst am: 01.12.2011, 14:23 Titel:

Hallo yiiit,

Nein, Deine DGL sieht nicht richtig aus. dot(3) und dot(4) sollten sich unterscheiden.

Gruß, Jan

yiiit

Themenstarter

Forum-Newbie


	Beiträge: 8

	Anmeldedatum: 21.10.11

	Wohnort: Osnabrück

	Version: ---

Verfasst am: 01.12.2011, 14:28 Titel:

ja, klar, logisch - damit y(4) die Geschwindigkeit in y-Richtung ergibt, muss am ende von dot(4) natürlich y(2) statt y(1) stehen.

edit: so besser? und was ist mit meinen anderen aussagen? Wink

danke!

yiiit

Themenstarter

Forum-Newbie


	Beiträge: 8

	Anmeldedatum: 21.10.11

	Wohnort: Osnabrück

	Version: ---

Verfasst am: 01.12.2011, 16:35 Titel:

Code:

function [] = main

[T,Y] = ode45('DGL',[0 10000], [10000 0 -8*cos(10) 8*sin(10)]);

plot(Y(:,1), Y(:,2))

axis([0 10000 -1 1])

Funktion ohne Link?

So sieht meine main aus.. warum erhalte ich, wenn ich die Zeile "axis(.." auskommentiere, eine Anzeige im Plot und wenn ich sie drin lasse, keine Ausgabe? noch dazu ist die Ausgabe im Plot ja nicht das, was man erwarten würde, da die Geschwindigkeit in y-Richtung ja positiv und nicht negativ ist..

Wenn ich cos(0) und sin(0) nehme, dann erhalte ich die zu erwartende gerade Linie, die für x-Werte von 0 bis 10000 einen y-Wert habe - aber nur mit auskommentierter "axis(.."-Zeile - sonst habe ich nur Wertepaare bis x=5000?

Was mache ich denn mit der tspan [0 10000]?

edit: irgendwas stimmt nicht.. meine Bahnkurve soll für größere Winkel zunächst positive y-Werte annehmen und mit der Zeit halt erst ansteigen und irgendwann wieder abfallen, sprich wie ein schräger Wurf aussehen - stattdessen sieht es aber so aus als würde die Raumkapsel abgestoßen werden, teilweise sind aber auch die y-Werte zunächst negativ, obwohl der Winkel ja schräg nach oben zeigt?

Jan S

Moderator


	Beiträge: 11.057

	Anmeldedatum: 08.07.10

	Wohnort: Heidelberg

	Version: 2009a, 2016b

Verfasst am: 01.12.2011, 20:50 Titel:

Hallo yiiit,

Einerseits könnte das ein Vorzeichen-Fehler sein, wenn z.B. die Richtung der Gravitation nach oben zeigen würde (eventuell g=-9.81?). Andererseits ist dies sehr ungewöhnlich:

Code:

... [10000 0 -8*cos(10) 8*sin(10)]);

Funktion ohne Link?

Das Argument von SIN und COS wird im Bogenmaß angegeben, also mit Werten zwischen 9 und 2*Pi. Meinst Du "cos(10 * pi / 180)" oder "cosd(10)"?

Gruß, Jan

yiiit

Themenstarter

Forum-Newbie


	Beiträge: 8

	Anmeldedatum: 21.10.11

	Wohnort: Osnabrück

	Version: ---

Verfasst am: 01.12.2011, 23:35 Titel:

danke jan!

die ergebnisse mit cosd und sind und mit -g sehen hervorragend aus!

Jan S

Moderator


	Beiträge: 11.057

	Anmeldedatum: 08.07.10

	Wohnort: Heidelberg

	Version: 2009a, 2016b

Verfasst am: 02.12.2011, 00:22 Titel:

Hallo yiiit,

Gerne!

Gruß, Jan

yiiit

Themenstarter

Forum-Newbie


	Beiträge: 8

	Anmeldedatum: 21.10.11

	Wohnort: Osnabrück

	Version: ---

Verfasst am: 07.12.2011, 14:45 Titel: Re: Lösung von DGL mit ode45

ich bin's nochmal!

jetzt soll auch noch die reibung berücksichtigt werden, die DGLs sehen jetzt so aus:

$\dot{v_{x}} = g \cdot \frac{r_{E}^{2}}{(x^{2}+y^{2})^{\frac{3}{2}}} \cdot x \cdot \frac{1}{2}A\rho_{L}\cdot v_{x}^{2}\cdot e^{-\frac{\sqrt{x^{2}+y^{2}}-R}{Hm}}$

$\dot{x} = v_{x}$

$\dot{v_{y}} = g \cdot \frac{r_{E}^{2}}{(x^{2}+y^{2})^{\frac{3}{2}}} \cdot y \cdot \frac{1}{2}A\rho_{L}\cdot v_{y}^{2}\cdot e^{-\frac{\sqrt{x^{2}+y^{2}}-R}{Hm}}$

$\dot{y} = v_{y}$

Dabei stellt ${\sqrt{x^{2}+y^{2}}-R$ die Höhe über dem Meeresspiegel dar, wenn ich mich nicht vertan habe..

Meine zugehörige Funktion dazu lautet:

Code:

function dot = DGL_e(t,y)

R = 6371;
g = 9.81/1000;

dot = [0;0;0;0];

dot(1) = y(3);
dot(2) = y(4);
dot(3) = -g.*R.^2./((y(1).^2+y(2).^2).^(3/2)).*y(1)+(5.*1.2.*y(3).^2.*exp(-(sqrt(y(1).^2+y(2).^2)-R)./8)./2000);
dot(4) = -g.*R.^2./((y(1).^2+y(2).^2).^(3/2)).*y(2)+(5.*1.2.*y(4).^2.*exp(-(sqrt(y(1).^2+y(2).^2)-R)./8)./2000);

Funktion ohne Link?

aber wenn ich das jetzt mit dieser Funktion auswerte:

Code:

function [] = main(a)

t = (0:pi/100:2*pi);
x = 6371 * sin(t);
y = 6371 * cos(t);

[T,Y] = ode45('DGL_e',[0 4000], [10000 0 -8.*cosd(a) 8.*sind(a)]);

Y(:,1:2)

plot(Y(:,1), Y(:,2))
hold on

axis([-10000 10000 -10000 10000])

plot(x,y)
% axis equal

Funktion ohne Link?

wo ich mir ja zur Hilfe die Werte ausgeben lassen, dann sieht man, dass die Werte "Not a Number" sind bei der Eingabe a = 5 zum Beispiel.. und ich erhalte kein vernünftiges Ergebnis..

dabei sollte man wieder ein Trajektorie sehen, bei der sich die Bahnkurve aufgrund der Reibung zu höheren Winkeln verschiebt als ohne, allerdings sieht man so höchstens, dass von der Erde abgeprallt wird.. ich brauche Hilfe!

edit: weiteres problem: wenn ich für a 0 einsetze, dann erhalte ich immer eine Warnung, auch bei der Berechnung ohne Reibung und das Ergebnis ist auch nicht das erwartete.. ab a=2 sind die Ergebnisse ohne Reibung perfekt.


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 - 2025 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.