Lösung gekoppelter Differentialgleichungen

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

Partner:

Forum

Option

[Erweitert]

• Diese Seite per Mail weiterempfehlen

Gehe zu:

paul101

Forum-Newbie


	Beiträge: 9

	Anmeldedatum: 03.01.10

	Wohnort: ---

	Version: ---

Verfasst am: 20.01.2010, 15:05 Titel: Lösung gekoppelter Differentialgleichungen

Hallo!

Ic habe diese Frage auch schon im octave-forum gestellt, aber da dort niemand geantwortet hat verschiebe ich sie mal hier hin.

Ich schreibe zur Zeit Facharbeit in Physik.
Ich bin gerade beim Berechnen eines Trebuchets
Dabei wird der Lagrange-Formalismus verwendet, ist aber nicht so wichtig...

Auf jeden Fall komme ich am Schluss auf diese drei Bewegungsgleichungen:

phi, sigma und theta sind jeweils Funktionen von der Zeit, die anderen Variablen sind kontstant.

jetzt zu meiner Frage:
gibt es eine Möglichkeit diese Gleichungen in Matlab/Octave zu lösen(ich bin ein ziemlicher Noob Confused

)?

Es wäre echt genial wenn mir jemand helfen würde Very Happy

Gruß,
Paul

Harald

Forum-Meister


	Beiträge: 24.499

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 20.01.2010, 15:19 Titel:

Hallo,

da hast du dir was vorgenommen.
Mehrere Schritte:
1. Schreibe das System in der Form M * dx/dt = f(x, t)
a) Löse die Gleichungen soweit möglich nach den zweiten Ableitungen auf, d.h. eine Seite enthält lineare Terme in den zweiten Ableitungen, die andere beliebig komplizierte Ausdrücke aus Funktionen und ersten Ableitungen.
b) setze t = theta', s = sigma' und d = delta' und stelle das Differentialgleichungssystem in den nun 6 Variablen auf.
c) Bringe das in die oben besagte Form.

2. Wende z.B. ode45 (in MATLAB) auf das Problem an. Wichtig: M muss als Parameter Mass über odeset übergeben werden.

Viel Erfolg,
Harald

paul101

Themenstarter

Forum-Newbie


	Beiträge: 9

	Anmeldedatum: 03.01.10

	Wohnort: ---

	Version: ---

Verfasst am: 20.01.2010, 20:09 Titel:

Danke für die Antwort!

Also bis 1a kapier ichs noch.
aber dei 1b setzt ich dann schon aus...
Was ist denn ein Differentialgleichungssystem?
Ist das sowas in der Art:

mit ph=phi', si=sigma' und th=theta'

wegen dem ode45:
das hier heißt auch ode45 und ist aber anscheinend für Octave:
kann ich das verwenden?

Gruß,
Paul

Harald

Forum-Meister


	Beiträge: 24.499

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 20.01.2010, 20:57 Titel:

Hallo,

du bist auf der richtigen Spur. Nur dass du 6 solche DGLen brauchst.
Differentialgleichungssystem = mehrere miteinander verknüpfte DGLen, wie unten von dir geschrieben.

Ein einfacheres Beispiel:
x'' = f(x', x, t)
als System erster Ordnung, mit u = x'
x' = u
u' = f(u, x, t)

Bei einer Facharbeit sollte dich auch der betreuende Lehrer da ein wenig unterstützen. Was du da hast, würde manch einen Studenten überfordern Wink

Grüße,
Harald

Harald

Forum-Meister


	Beiträge: 24.499

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 21.01.2010, 12:18 Titel:

Hallo,

noch ein Hinweis ... hier dürften nützliche Informationen zum Thema zu finden sein:
http://www.mathworks.com/support/tech-notes/1500/1510.html

Grüße,
Harald

paul101

Themenstarter

Forum-Newbie


	Beiträge: 9

	Anmeldedatum: 03.01.10

	Wohnort: ---

	Version: ---

Verfasst am: 21.01.2010, 17:28 Titel:

hallo,
ich habe heute mit meinem Physiklehrer darüber geredet.
So wie ich das jetzt verstanden habe muss ich folgendes tun:

1. Die 3 Gleichungen nach phidotdot, sigmadotdot und thetadotdot auflösen
2. Diese zweiten Ableitungen jeweils aus den 3 Gleichungen so rauseliminieren, dass ich dann 3 Gleichungen in der Form

phidotdot=f1(phidot, phi, sigmadot, sigma, thetadot, theta)
sigmadotdot=f2(phidot, phi, sigmadot, sigma, thetadot, theta)
thetadotdot=f3(phidot, phi, sigmadot, sigma, thetadot, theta)
habe
3. Das ganze so ähnlich als funktion schreiben:

Code:

function dy = f(y,t)
dy = zeros(6,1);

dy(1) = y(2);
dy(2) = f1(y1, y2, y3, y4, y5, y6);
dy(3) = y(4);
dy(4) = f2(y1, y2, y3, y4, y5, y6);
dy(5) = y(6);
dy(6) = f3(y1, y2, y3, y4, y5, y6);
enfunction

Funktion ohne Link?

4. Das dann an den Integrator übergeben

Gruß,
Paul

Harald

Forum-Meister


	Beiträge: 24.499

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 21.01.2010, 17:45 Titel:

Hallo,

Genau das ist der Weg.
2. Dürfte aber etwas schwierig werden, da deine Gleichungen recht komplex und schwer nach den zweiten Ableitungen aufzulösen sind. Hast du die Symbolic Math Toolbox oder äquivalent zur Verfügung?

Edit: Falls du die Gleichungen selber hergeleitet hast, solltest du sie auf jeden Fall nochmal überprüfen. Wäre ja ärgerlich, wenn du das alles nochmal durchgehen müsstest.

Grüße,
Harald

paul101

Themenstarter

Forum-Newbie


	Beiträge: 9

	Anmeldedatum: 03.01.10

	Wohnort: ---

	Version: ---

Verfasst am: 21.01.2010, 19:03 Titel:

Hallo,

ich habe die Gleichungen zwar selber, aber nicht von Hand hergeleitet. Dazu habe ich Derive in der Version 4(Die hab ich von der Schule) benutzt. Das war zwar anstrengend aber sicher weniger anstrengend als die von Hand herzuleiten...aber ich werde es auf jeden Fall noch einmal kontrollieren.

Ich habe leider nur Octave, deswegen bringt mir diese Toolbox glaub ich nichts...
Weißt du ob es etwas vergleichbares für Octave gibt?

Gruß,
Paul

paul101

Themenstarter

Forum-Newbie


	Beiträge: 9

	Anmeldedatum: 03.01.10

	Wohnort: ---

	Version: ---

Verfasst am: 24.01.2010, 16:34 Titel:

Hallo

Also erstmal vielen, vielen Dank für die tolle Hilfe!
Ich habe die Gleichungen jetzt mit Maxima aufgelöst, allerdings stimmt das vom Skript ausgespuckte Ergebnis nicht so wirklich...
Das Skipt soll die Bewegung des einzelnen Punkte, also Projektil etc, im Koordinatensystem anzeigen.
Bis zu einem bestimmten Zeitpunkt stimmt die Bewegung so einigermaßen, bis auf die (Im Plot blaue) Bewegung des Gegengewichts. Aber dann bricht die numerische Berechnung irgendwie ab weil alles zu schnell geht, so verstehe ich das aufjedenfall...
So schnell sollte sich das aber eigentlich nicht bewegen, deswegen frage ich mich jetzt ob das am Code liegt oder an den Gleichungen.
Könnt ihr mir sagen ob der Code falsch ist?

Code:

clear

printf("\n")

%konstanten:
global g=9.81;
global l1=input("l1 in m: ");
global l2=input("l2 in m: ");
global l3=input("l3 in m: ");
global l4=input("l4 in m: ");
global m1=input("m1 in kg: ");
global m2=input("m2 in kg: ");

%differentialgleichung
function dy= f(t,y)

global g
global l1
global l2
global l3
global l4
global m1
global m2

dy = zeros(6,1);

%ddotphi:
dy(1) = y(2);
dy(2) = -((y(2)^2*l1^2*m1*cos(y(1))-g*l1*m1)*sin(y(1))*sin(y(1)-y(3))^2+((-y(2)^2*l1^2*m1*sin(y(1))^2+y(2)^2*l1^2*m1*cos(y(1))^2-g*l1*m1*cos(y(1)))*cos(y(1)-y(3))+(-y(4)^2+2*y(2)*y(4)-y(2)^2)*l1*l3*m1*cos(y(1)))*sin(y(1)-y(3))-y(2)^2*l1^2*m1*cos(y(1))*sin(y(1))*cos(y(1)-y(3))^2+(y(4)^2-2*y(2)*y(4)+y(2)^2)*l1*l3*m1*sin(y(1))*cos(y(1)-y(3))+(y(2)^2*l2^2*m2*cos(y(1))+g*l2*m2)*sin(y(1))*sin(y(5)+y(1))^2+((y(2)^2*l2^2*m2*sin(y(1))^2+y(2)^2*l2^2*m2*cos(y(1))^2+g*l2*m2*cos(y(1)))*cos(y(5)+y(1))+(-y(6)^2-2*y(2)*y(6)-y(2)^2)*l2*l4*m2*cos(y(1)))*sin(y(5)+y(1))+y(2)^2*l2^2*m2*cos(y(1))*sin(y(1))*cos(y(5)+y(1))^2+(y(6)^2+2*y(2)*y(6)-y(2)^2)*l2*l4*m2*sin(y(1))*cos(y(5)+y(1))+(g*l1*m1-g*l2*m2)*sin(y(1)))/(l1^2*m1*sin(y(1))^2*sin(y(1)-y(3))^2+2*l1^2*m1*cos(y(1))*sin(y(1))*cos(y(1)-y(3))*sin(y(1)-y(3))+l1^2*m1*cos(y(1))^2*cos(y(1)-y(3))^2+l2^2*m2*sin(y(1))^2*sin(y(5)+y(1))^2+2*l2^2*m2*cos(y(1))*sin(y(1))*cos(y(5)+y(1))*sin(y(5)+y(1))+l2^2*m2*cos(y(1))^2*cos(y(5)+y(1))^2-l2^2*m2-l1^2*m1) ;

%ddottheta:
dy(3) = y(4);
dy(4) = -(((y(4)^2-2*y(2)*y(4)+2*y(2)^2)*l1^2*l3*m1*cos(y(1))-g*l1*l3*m1)*sin(y(1))*sin(y(1)-y(3))^2+(((-y(4)^2+2*y(2)*y(4)-2*y(2)^2)*l1^2*l3*m1*sin(y(1))^2+(y(4)^2-2*y(2)*y(4)+2*y(2)^2)*l1^2*l3*m1*cos(y(1))^2-g*l1*l3*m1*cos(y(1)))*cos(y(1)-y(3))+(-g*l2^2-g*l1*l2)*m2*sin(y(1))^2*sin(y(5)+y(1))^2+(((y(2)^2*l1*l2^2*m2*cos(y(1))^2+(-2*g*l2^2-g*l1*l2)*m2*cos(y(1)))*sin(y(1))-y(2)^2*l1*l2^2*m2*sin(y(1))^3)*cos(y(5)+y(1))+(y(6)^2+2*y(2)*y(6)+y(2)^2)*l1*l2*l4*m2*cos(y(1))*sin(y(1)))*sin(y(5)+y(1))+(-y(2)^2*l1*l2^2*m2*cos(y(1))*sin(y(1))^2+y(2)^2*l1*l2^2*m2*cos(y(1))^3-g*l2^2*m2*cos(y(1))^2)*cos(y(5)+y(1))^2+(-y(6)^2-2*y(2)*y(6)+y(2)^2)*l1*l2*l4*m2*sin(y(1))^2*cos(y(5)+y(1))+(g*l1*l2*m2-g*l1^2*m1)*sin(y(1))^2+(((-y(4)^2+2*y(2)*y(4)-y(2)^2)*l1*l3^2-y(2)^2*l1^3)*m1-y(2)^2*l1*l2^2*m2)*cos(y(1))+g*l2^2*m2+g*l1^2*m1)*sin(y(1)-y(3))+(-y(4)^2+2*y(2)*y(4)-2*y(2)^2)*l1^2*l3*m1*cos(y(1))*sin(y(1))*cos(y(1)-y(3))^2+(((-y(2)^2*l1*l2^2*m2*cos(y(1))^2-g*l1*l2*m2*cos(y(1)))*sin(y(1))-y(2)^2*l1*l2^2*m2*sin(y(1))^3)*sin(y(5)+y(1))^2+((-3*y(2)^2*l1*l2^2*m2*cos(y(1))*sin(y(1))^2-y(2)^2*l1*l2^2*m2*cos(y(1))^3-g*l1*l2*m2*cos(y(1))^2)*cos(y(5)+y(1))+(y(6)^2+2*y(2)*y(6)+y(2)^2)*l1*l2*l4*m2*cos(y(1))^2)*sin(y(5)+y(1))-2*y(2)^2*l1*l2^2*m2*cos(y(1))^2*sin(y(1))*cos(y(5)+y(1))^2+(-y(6)^2-2*y(2)*y(6)+y(2)^2)*l1*l2*l4*m2*cos(y(1))*sin(y(1))*cos(y(5)+y(1))+((g*l1*l2*m2-g*l1^2*m1)*cos(y(1))+y(2)^2*l1*l2^2*m2+((y(4)^2-2*y(2)*y(4)+y(2)^2)*l1*l3^2+y(2)^2*l1^3)*m1)*sin(y(1)))*cos(y(1)-y(3))+(y(2)^2*l2^2*l3*m2*cos(y(1))+g*l2*l3*m2)*sin(y(1))*sin(y(5)+y(1))^2+((y(2)^2*l2^2*l3*m2*sin(y(1))^2+y(2)^2*l2^2*l3*m2*cos(y(1))^2+g*l2*l3*m2*cos(y(1)))*cos(y(5)+y(1))+(-y(6)^2-2*y(2)*y(6)-y(2)^2)*l2*l3*l4*m2*cos(y(1)))*sin(y(5)+y(1))+y(2)^2*l2^2*l3*m2*cos(y(1))*sin(y(1))*cos(y(5)+y(1))^2+(y(6)^2+2*y(2)*y(6)-y(2)^2)*l2*l3*l4*m2*sin(y(1))*cos(y(5)+y(1))+(g*l1*l3*m1-g*l2*l3*m2)*sin(y(1)))/(l1^2*l3*m1*sin(y(1))^2*sin(y(1)-y(3))^2+2*l1^2*l3*m1*cos(y(1))*sin(y(1))*cos(y(1)-y(3))*sin(y(1)-y(3))+l1^2*l3*m1*cos(y(1))^2*cos(y(1)-y(3))^2+l2^2*l3*m2*sin(y(1))^2*sin(y(5)+y(1))^2+2*l2^2*l3*m2*cos(y(1))*sin(y(1))*cos(y(5)+y(1))*sin(y(5)+y(1))+l2^2*l3*m2*cos(y(1))^2*cos(y(5)+y(1))^2-l2^2*l3*m2-l1^2*l3*m1) ;

%ddotsigma:
dy(5) = y(6);
dy(6) = (((g*l1*l2+g*l1^2)*m1*sin(y(1))^2*sin(y(5)+y(1))+(y(2)^2*l1^2*l2*m1*sin(y(1))^3+(g*l1*l2*m1*cos(y(1))-y(2)^2*l1^2*l2*m1*cos(y(1))^2)*sin(y(1)))*cos(y(5)+y(1))+(y(2)^2*l1^2*l4*m1*cos(y(1))-g*l1*l4*m1)*sin(y(1)))*sin(y(1)-y(3))^2+(((y(2)^2*l1^2*l2*m1*sin(y(1))^3+(y(2)^2*l1^2*l2*m1*cos(y(1))^2+(g*l1*l2+2*g*l1^2)*m1*cos(y(1)))*sin(y(1)))*sin(y(5)+y(1))+(3*y(2)^2*l1^2*l2*m1*cos(y(1))*sin(y(1))^2-y(2)^2*l1^2*l2*m1*cos(y(1))^3+g*l1*l2*m1*cos(y(1))^2)*cos(y(5)+y(1))-y(2)^2*l1^2*l4*m1*sin(y(1))^2+y(2)^2*l1^2*l4*m1*cos(y(1))^2-g*l1*l4*m1*cos(y(1)))*cos(y(1)-y(3))+(y(4)^2-2*y(2)*y(4)+y(2)^2)*l1*l2*l3*m1*cos(y(1))*sin(y(1))*sin(y(5)+y(1))+(y(4)^2-2*y(2)*y(4)+y(2)^2)*l1*l2*l3*m1*cos(y(1))^2*cos(y(5)+y(1))+(-y(4)^2+2*y(2)*y(4)-y(2)^2)*l1*l3*l4*m1*cos(y(1)))*sin(y(1)-y(3))+((y(2)^2*l1^2*l2*m1*cos(y(1))*sin(y(1))^2+y(2)^2*l1^2*l2*m1*cos(y(1))^3+g*l1^2*m1*cos(y(1))^2)*sin(y(5)+y(1))+2*y(2)^2*l1^2*l2*m1*cos(y(1))^2*sin(y(1))*cos(y(5)+y(1))-y(2)^2*l1^2*l4*m1*cos(y(1))*sin(y(1)))*cos(y(1)-y(3))^2+((-y(4)^2+2*y(2)*y(4)-y(2)^2)*l1*l2*l3*m1*sin(y(1))^2*sin(y(5)+y(1))+(-y(4)^2+2*y(2)*y(4)-y(2)^2)*l1*l2*l3*m1*cos(y(1))*sin(y(1))*cos(y(5)+y(1))+(y(4)^2-2*y(2)*y(4)+y(2)^2)*l1*l3*l4*m1*sin(y(1)))*cos(y(1)-y(3))+((y(6)^2+2*y(2)*y(6)+2*y(2)^2)*l2^2*l4*m2*cos(y(1))+g*l2*l4*m2)*sin(y(1))*sin(y(5)+y(1))^2+(((-y(6)^2-2*y(2)*y(6)+2*y(2)^2)*l2^2*l4*m2*sin(y(1))^2+(y(6)^2+2*y(2)*y(6)+2*y(2)^2)*l2^2*l4*m2*cos(y(1))^2+g*l2*l4*m2*cos(y(1)))*cos(y(5)+y(1))+(g*l2^2*m2-g*l1*l2*m1)*sin(y(1))^2+(((-y(6)^2-2*y(2)*y(6)-y(2)^2)*l2*l4^2-y(2)^2*l2^3)*m2-y(2)^2*l1^2*l2*m1)*cos(y(1))-g*l2^2*m2-g*l1^2*m1)*sin(y(5)+y(1))+(-y(6)^2-2*y(2)*y(6)+2*y(2)^2)*l2^2*l4*m2*cos(y(1))*sin(y(1))*cos(y(5)+y(1))^2+((g*l2^2*m2-g*l1*l2*m1)*cos(y(1))+((y(6)^2+2*y(2)*y(6)-y(2)^2)*l2*l4^2-y(2)^2*l2^3)*m2-y(2)^2*l1^2*l2*m1)*sin(y(1))*cos(y(5)+y(1))+(g*l1*l4*m1-g*l2*l4*m2)*sin(y(1)))/(l1^2*l4*m1*sin(y(1))^2*sin(y(1)-y(3))^2+2*l1^2*l4*m1*cos(y(1))*sin(y(1))*cos(y(1)-y(3))*sin(y(1)-y(3))+l1^2*l4*m1*cos(y(1))^2*cos(y(1)-y(3))^2+l2^2*l4*m2*sin(y(1))^2*sin(y(5)+y(1))^2+2*l2^2*l4*m2*cos(y(1))*sin(y(1))*cos(y(5)+y(1))*sin(y(5)+y(1))+l2^2*l4*m2*cos(y(1))^2*cos(y(5)+y(1))^2-l2^2*l4*m2-l1^2*l4*m1) ;
endfunction

%Anfangswinkel phi0:
phiadeg=input("phianfang in Grad: ");
global phi0=phiadeg/180*pi;

%Anfangswinkel theta0:
global theta0=phi0-pi/2;

%Anfangsiwnkel sigma0:
global sigma0=pi - phi0;

%Anfangsbedingungen:
anf =[phi0, 0, theta0 ,0 ,sigma0, 0];

%Lösen der Differentialgleichung:
[t,L]=ode45("f",[0 5],anf);

%position des Projektils:
x1=-l1*sin(L(:,1))+l3*sin(L(:,1)-L(:,3));
y1=l1*cos(L(:,1))-l3*cos(L(:,3)-L(:,1));

%Position des GG:
x2=l2*sin(L(:,1))-l4*sin(L(:,1)+L(:,5));
y2=-l2*cos(L(:,1))+l4*cos(L(:,1)+L(:,5));

%Position des vorderen Ende des Arms:
x3=-l1*sin(L(:,1));
y3=l1*cos(L(:,1));

%Position des hinteren Ende des Arms:
x4=l2*sin(L(:,1));
y4=-l2*cos(L(:,1));

%Plot:
figure
subplot(2,1,1);
plot(t,L(:,1),"r-");
title("phi");
grid on

subplot(2,1,2);
plot(t,L(:,2),"r-");
title("phidot");
grid on

figure;
plot(x1,y1,"-r", x2,y2,"-b", x3,y3,"-m", x4,y4,"-c", 0,0,"og;;");
title("???");
xlabel("x in m");
ylabel("y in m");
grid on

Funktion ohne Link?

Gruß,
Paul

Harald

Forum-Meister


	Beiträge: 24.499

	Anmeldedatum: 26.03.09

	Wohnort: Nähe München

	Version: ab 2017b

Verfasst am: 24.01.2010, 16:55 Titel:

Hallo,

ist natürlich jetzt schwierig, das zu sagen. Ob irgendwo in den Gleichungen ein Fehler drin ist, ist schwer bis unmöglich zu sagen. Wenn das Verhalten bis zu einem gewissen Punkt stimmt, liegen die Probleme aber vermutlich woanders.

Vorschläge, bezogen auf MATLAB:
- Toleranz heruntersetzen (mit odeset)
- Testen anderer Löser (z.B. ode23s)

Zitat:

Bis zu einem bestimmten Zeitpunkt stimmt die Bewegung so einigermaßen, bis auf die (Im Plot blaue) Bewegung des Gegengewichts.

Dann rechne am besten diese Gleichung nochmal nach. In solchen Fällen ist es oft einfacher, es tatsächlich nochmal durchzurechnen und dann zu vergleichen als Fehler zu suchen.

Grüße,
Harald


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.