| |
| Autor |
  Einmasseschwinger DGL lösen und plotten |
|
Tim80
Ehemals Aktiv 
Dabei seit: 11.05.2015
Mitteilungen: 47
 |  Themenstart: 2017-05-24
|
Hallo,
ich möchte an einem einfachen Beispiel testen wie man in Matlab Differentialgleichungen lösen und Ergebnisse plotten kann.
Dazu habe ich mir einen einfachen Einmasseschwinger, oder auch Feder-Masse-Schwinger genannt, ausgesucht.
m*q^**+d*q^*+c*q=F
q^**=F/m-(d/m)*q^*-(c/m)*q
Die DGL steht nun fest und wurde nach der Beschleunigung q^** umgestellt.
Ich möchte im Folgenden den Weg, die Geschwindigkeit und die Beschleunigung der Masse über die Zeit plotten.
Das Plotten des Wegs und der Geschwindigkeit hat geklappt, nur bei der Beschleunigung funktioniert es nicht.
Hier wird der Wert irgendwie nicht aktualisiert. Ich habe wirklich keine Ahnung woran das liegt.
Es wäre super wenn mir jemand weiterhelfen könnte.
Grüße
Tim
clear all, close all, clc, clf;
global m c d g F;
m = 1; %Masse
c = 0; %Steifigkeitskoeffizient
d = 0; %Dämpfungskoeffizient
g = 9.81; %Gravitationskonstante
F = m*g; %Gewichtskraft
%Anfangsbedingungen
q0=0;
v0=0;
a0=0;
AB=[q0,v0,a0];
tspan = [0:0.01:10]; %Zeitvektor
[t,q] = ode45(@FCN_EMS,tspan,AB);
q(:,3)=(F/m)-(d/m)*q(2)-(c/m)*q(1);
%F=50*abs(sin(2*t*pi));
%F=1*9.81;
%q(:,3)=(F/1)-(0.5/1)*q(2)-(100/1)*q(1);
figure(1);
subplot(311)
Plot_q = plot(t,q(:,1)); %Plot der Weg über Zeit
hold on,grid on
xlabel('Zeit in Sekunden') %X-Achsenbeschriftung
ylabel('Weg in m') %Y-Achsenbeschriftung
subplot(312)
Plot_v = plot(t,q(:,2)); %Plot der Geschwindigkeit über Zeit
hold on,grid on
xlabel('Zeit in Sekunden') %X-Achsenbeschriftung
ylabel('Geschwindigkeit in m/s') %Y-Achsenbeschriftung
subplot(313)
Plot_a = plot(t,q(:,3)); %Plot der Beschleunigung über Zeit
hold on,grid on
xlabel('Zeit in Sekunden') %X-Achsenbeschriftung
ylabel('Beschleunigung in m/s^2') %Y-Achsenbeschriftung
function [ dq ] = FCN_EMS( t,q )
global m c d g F;
%q(1) = AB1 --> q0 Höhe zum Zeitpunkt t=0s
%q(2) = AB2 --> v0 Geschwindigkeit zum Zeitpunkt t=0s
%q(3) = AB3 --> a0 Beschleunigung zum eitpunkt t=0s
dq = zeros(3,1); %Vordefinition Zustandsvektor
dq(1) = q(2); %Geschwindigkeit
dq(2) = (F/m)-(d/m)*q(2)-(c/m)*q(1); %Beschleunigung
end
|
 Profil
|
Tim80
Ehemals Aktiv 
Dabei seit: 11.05.2015
Mitteilungen: 47
| Beitrag No.1, vom Themenstarter, eingetragen 2017-05-24
|
Das Problem mit dem Plotten der Beschleunigung habe ich lösen können.
Es ist zwar nicht sonderlich elegant, aber dennoch funktioniert es.
Diese zusätzliche Zeile habe ich außerhalb der Funktion hinzugefügt:
q(:,3)=q(:,2)./t;
Ich bin aber natürlich für bessere Lösung offen :-) (Unten der neue Quellcode)
Jetzt habe ich aber noch ein anderes Problem:
Ich möchte anstatt der Gravitationskraft eine zeitabhängige schwellende Kraft wirken lassen (sinusförmig).
Diese habe ich mir wie folgt vorgestellt: F=10*sin(t.*pi); wobei t die Zeit ist.
Wie kann ich eine solche zeitabhängige Kraft ins System einbringen? Wo muss ich diese definieren?
Ich würde mich über eure Hilfe sehr freuen!
Grüße
PS. Wie kann man eigentlich den Matlab-Quellcode "schön" in eine Forenbeitrag kopieren/einfügen? In anderen Threads habe ich das schon gesehen, hab aber nicht rausfinden können wie es geht...
clear all, close all, clc, clf;
global m c d g F;
m = 2; %Masse
c = 0; %Steifigkeitskoeffizient
d = 0; %Dämpfungskoeffizient
g = 9.81; %Gravitationskonstante
F = m*g; %Externe Kraft
%Anfangsbedingungen
q0=0;
v0=0;
a0=0;
AB=[q0,v0,a0];
tspan = transp([0:0.01:10]); %Zeitvektor
[t,q] = ode45(FCN_EMS,tspan,AB);
q(:,3)=q(:,2)./t;
figure(1);
subplot(311)
Plot_q = plot(t,q(:,1)); %Plot der Weg über Zeit
hold on,grid on
xlabel('Zeit in Sekunden') %X-Achsenbeschriftung
ylabel('Weg in m') %Y-Achsenbeschriftung
subplot(312)
Plot_v = plot(t,q(:,2)); %Plot der Geschwindigkeit über Zeit
hold on,grid on
xlabel('Zeit in Sekunden') %X-Achsenbeschriftung
ylabel('Geschwindigkeit in m/s') %Y-Achsenbeschriftung
subplot(313)
Plot_a = plot(t,q(:,3)); %Plot der Beschleunigung über Zeit
hold on,grid on
xlabel('Zeit in Sekunden') %X-Achsenbeschriftung
ylabel('Beschleunigung in m/s^2') %Y-Achsenbeschriftung
function [ dq ] = FCN_EMS( t,q )
global m c d F;
%q(1) = AB1 --> q0 Höhe zum Zeitpunkt t=0s
%q(2) = AB2 --> v0 Geschwindigkeit zum Zeitpunkt t=0s
%q(3) = AB3 --> a0 Beschleunigung zum Zeitpunkt t=0s
dq = zeros(3,1); %Vordefinition Zustandsvektor
dq(1) = q(2); %Geschwindigkeit
dq(2) = (F/m)-(d/m)*q(2)-(c/m)*q(1); %Beschleunigung
end
|
Profil
|
majoka
Senior 
Dabei seit: 25.02.2014
Mitteilungen: 810
| Beitrag No.2, eingetragen 2017-05-24
|
Die Beschleunigung sollte sich doch einfach so ergeben
\sourceon Matlab
q(:,3) = (F(t)/m)-(d/m)*q(:,2)-(c/m)*q(:,1);
\sourceoff
Die zeitliche Änderung von F kannst Du ganz einfach mit in die Funktion aufnehmen.
\sourceon Matlab
%Anfangsbedingungen
q0=0;
v0=0;
%a0=0;
AB=[q0,v0];
tspan = [0 20]; %Zeitvektor
F = @(t) 10*sin(t*pi);
FCN_EMS = @(t,q) [q(2);(F(t)/m)-(d/m)*q(2)-(c/m)*q(1)];
[t,q] = ode45(FCN_EMS,tspan,AB);
q(:,3) = (F(t)/m)-(d/m)*q(:,2)-(c/m)*q(:,1);
\sourceoff
Für die Einbindung des Codes drücke auf "Quelltextbereich" und ersetze dann "nameDerSprache" durch "Matlab".
|
Profil
|
Tim80
Ehemals Aktiv
Dabei seit: 11.05.2015
Mitteilungen: 47
| Beitrag No.3, vom Themenstarter, eingetragen 2017-05-29
|
\
Hallo majoka,
vielen Dank für deine Hilfe. Ich habe das Ganze in mein Programm integriert und jetzt funktioniert es auch für eine zeitlich veränderliche externe Kraft. Mein neuer Quellcode habe ich unten eingefügt.
Wenn nun aber beispielsweise nur die Gravitationskraft F_G = m*g wirken soll, funktioniert es leider nicht mehr, da diese nicht zeitabhängig ist.
Die Kraft darf dann eigentlich nicht mehr in der Funktion stehen und überall im Programm muss die Zeitabhängigkeit (t) bei der Kraft F(t) entfernt werden. Das ist natürlich aufwändig und auch nicht gewünscht.
Gibt es hierfür eine "schöne" Lösung?
Grüße Tim
\sourceon Matlab
clear all, close all, clc, clf;
global m c d g F; %Definition globale Variablen
m = 2; %Masse
c = 0; %Steifigkeitskoeffizient
d = 0; %Dämpfungskoeffizient
g = 9.81; %Gravitationskonstante
%Anfangsbedingungen
q0=0;
v0=0;
a0=0;
AB=[q0,v0,a0];
tspan = transp([0:0.1:20]); %Zeitvektor
[t,q] = ode45(@FCN_EMS,tspan,AB); %ODE45 Solver
q(:,3)=(F(t)/m)-(d/m)*q(:,2)-(c/m)*q(:,1); %Beschleunigung (für Plot)
figure(1);
subplot(311) %Plot der Weg über Zeit
Plot_q = plot(t,q(:,1));
hold on,grid on
xlabel('Zeit in Sekunden') %X-Achsenbeschriftung
ylabel('Weg in m') %Y-Achsenbeschriftung
subplot(312) %Plot der Geschwindigkeit über Zeit
Plot_v = plot(t,q(:,2));
hold on,grid on
xlabel('Zeit in Sekunden') %X-Achsenbeschriftung
ylabel('Geschwindigkeit in m/s') %Y-Achsenbeschriftung
subplot(313) %Plot der Beschleunigung über Zeit
Plot_a = plot(t,q(:,3));
hold on,grid on
xlabel('Zeit in Sekunden') %X-Achsenbeschriftung
ylabel('Beschleunigung in m/s^2') %Y-Achsenbeschriftung
figure(2)
Plot_F = plot(t,F(t)); %Plot der externen Kraft über Zeit
hold on,grid on
xlabel('Zeit in Sekunden') %X-Achsenbeschriftung
ylabel('Externe Kraft in N') %Y-Achsenbeschriftung
function [ dq ] = FCN_EMS( t,q )
global m c d F; %Definition/Aufruf globale Variablen
%q(1) = AB1 --> q0 Höhe zum Zeitpunkt t=0s
%q(2) = AB2 --> v0 Geschwindigkeit zum Zeitpunkt t=0s
%q(3) = AB3 --> a0 Beschleunigung zum Zeitpunkt t=0s
F = @(t) 10*sin(t*pi); %Externe Kraft
dq = zeros(3,1); %Vordefinition Zustandsvektor
dq(1) = q(2); %Geschwindigkeit
dq(2) = (F(t)/m)-(d/m)*q(2)-(c/m)*q(1); %Beschleunigung
end
\sourceoff
|
Profil
|
majoka
Senior
Dabei seit: 25.02.2014
Mitteilungen: 810
| Beitrag No.4, eingetragen 2017-05-29
|
Die Fehlermeldung sollte sich mit
\sourceon Matlab
g = 9.81;
F = @(t) m*g+0*t;
\sourceoff
vermeiden lassen.
|
Profil
|
Tim80
Ehemals Aktiv
Dabei seit: 11.05.2015
Mitteilungen: 47
| Beitrag No.5, vom Themenstarter, eingetragen 2017-05-30
|
Super, vielen Dank, damit hat es funktioniert!
Jetzt habe ich noch eine weitere Frage :
Ich möchte für die Beschleunigung die Anfangsbedingung a0 vorgeben, genauso wie ich es auch für q0 und v0 gemacht habe.
Ich habe versucht den Code zu erweitern, leider sind nun alle Werte für Weg und Geschwindigkeit Null. Ich weiß nicht wo der Fehler genau liegt, ich habe gedacht, dass ich nur den Vektor für die Anfangsbedingungen und die Funktion an sich erweitern muss.
Könntest du mir bitte dabei nochmal helfen?
LG
\sourceon Matlab
clear all, close all, clc, clf;
global m c d g F; %Definition globale Variablen
m = 2; %Masse
c = 0; %Steifigkeitskoeffizient
d = 0; %Dämpfungskoeffizient
g = 9.81; %Gravitationskonstante
%Anfangsbedingungen
q0=0;
v0=0;
a0=0;
AB=[q0,v0,a0];
tspan = [0:0.1:20]; %Zeitvektor
F = @(t) 10*sin(t*pi); %Externe Kraft
%F = @(t) m*g+0*t; %Gravitationskraft
FCN_EMS = @(t,q) [q(1);q(2);(F(t)/m)-(d/m)*q(2)-(c/m)*q(1)];
[t,q] = ode45(FCN_EMS,tspan,AB);
%q(:,3) = (F(t)/m)-(d/m)*q(:,2)-(c/m)*q(:,1);
figure(1);
subplot(311)
Plot_q = plot(t,q(:,1)); %Plot der Weg über Zeit
hold on,grid on
xlabel('Zeit in Sekunden') %X-Achsenbeschriftung
ylabel('Weg in m') %Y-Achsenbeschriftung
subplot(312)
Plot_v = plot(t,q(:,2)); %Plot der Geschwindigkeit über Zeit
hold on,grid on
xlabel('Zeit in Sekunden') %X-Achsenbeschriftung
ylabel('Geschwindigkeit in m/s') %Y-Achsenbeschriftung
subplot(313)
Plot_a = plot(t,q(:,3)); %Plot der Beschleunigung über Zeit
hold on,grid on
xlabel('Zeit in Sekunden') %X-Achsenbeschriftung
ylabel('Beschleunigung in m/s^2') %Y-Achsenbeschriftung
figure(2)
plotF = plot(t,F(t));
hold on,grid on
xlabel('Zeit in Sekunden') %X-Achsenbeschriftung
ylabel('Externe Kraft in N') %Y-Achsenbeschriftung
\sourceoff
|
Profil
|
majoka
Senior
Dabei seit: 25.02.2014
Mitteilungen: 810
| Beitrag No.6, eingetragen 2017-05-30
|
Ich habe mir es jetzt nicht im Detail angeschaut, aber bei einer DGL 2. Ordnung hat man doch nur zwei Anfangsbedingungen!?
Der Wert $a(0)$ ergibt sich doch aus der DGL.
|
Profil
|
Tim80
Ehemals Aktiv
Dabei seit: 11.05.2015
Mitteilungen: 47
| Beitrag No.7, vom Themenstarter, eingetragen 2017-05-30
|
Stimmt, Fehler meinerseits, da hab ich gepennt...
Danke für deine Hilfe!
LG
|
Profil
|
Tim80 hat die Antworten auf ihre/seine Frage gesehen.
Tim80 hat selbst das Ok-Häkchen gesetzt. |
|
All logos and trademarks in this site are property of their respective owner. The comments are property of their posters, all the rest © 2001-2023 by Matroids Matheplanet
This web site was originally made with PHP-Nuke, a former web portal system written in PHP that seems no longer to be maintained nor supported. PHP-Nuke is Free Software released under the GNU/GPL license.
Ich distanziere mich von rechtswidrigen oder anstößigen Inhalten, die sich trotz aufmerksamer Prüfung hinter hier verwendeten Links verbergen mögen.
Lesen Sie die
Nutzungsbedingungen,
die Distanzierung,
die Datenschutzerklärung und das Impressum.
[Seitenanfang]
|
Home / Seite ohne Frame
Wie unterstütze ich den Matheplaneten mit Geld?
