Nichtlineare Schwingungsdynamik¶
Übung - Aufgabe 1 ¶
Bestimmen Sie den Phasenplan eines Ein-Massen-Schwingers mit der Rückstellfunktion $$ f(x) = \sin(x) $$ analytisch.
Stellen Sie zudem die Näherungen der 1. Ordnung (Linearisierung), 3. und 5. Ordnung den Ergebnissen des Originalmodells mit verschiedenen Anfangsbedingungen gegenüber.
$$
J_A \ddot{\varphi} + mga \sin(\varphi) = 0
$$
In [1]:
% init octave
close all
clear
clc
set(0, "defaultlinelinewidth", 5.0); % Linienbreite fuer alle plots setzen
set(0, "defaultaxesfontsize", 16)
set(0, "defaulttextfontsize", 16)
Systemparameter
$$ J = 1.0,\qquad mga = 1.0,\qquad E_0 = 0.9 $$In [2]:
J = 1.0;
mga = 1.0;
E0 = 0.9;
Analytische Lösungen
$E_{pot}$
In [3]:
phi_max = acos(1-E0/mga);
phi = linspace(-phi_max,phi_max,100);
Epot_ex = mga*(1-cos(phi));
Epot_1O = 0.5*phi.^2;
Epot_3O = 0.5*phi.^2-1/24*phi.^4;
Epot_5O = 0.5*phi.^2-1/24*phi.^4+1/720*phi.^6;
$E_{pot}$ darstellen:
In [4]:
figure('position',[0 0 1000 500])
plot(phi,Epot_ex,'k','DisplayName','sin(\varphi)');
hold on;
plot(phi,Epot_1O,'b','DisplayName','\varphi');
plot(phi,Epot_3O,'g','DisplayName','\varphi - 1/6 \varphi^3');
plot(phi,Epot_5O,'r','DisplayName','\varphi - 1/6 \varphi^3 + 1/120 \varphi^5');
title('analytische Loesung')
ylabel('Epot')
xlabel('\varphi')
h = legend('location', 'eastoutside');
set(h, "interpreter", "tex");
Phasenplan
In [5]:
phid_ex_pos = + sqrt(2/J * (E0+mga*(cos(phi)-1)));
phid_ex_neg = - sqrt(2/J * (E0+mga*(cos(phi)-1)));
phid_1O_pos = + sqrt(2/J * (E0-mga*(0.5*phi.^2)));
phid_1O_neg = - sqrt(2/J * (E0-mga*(0.5*phi.^2)));
phid_3O_pos = + sqrt(2/J * (E0-mga*(0.5*phi.^2-1/24*phi.^4)));
phid_3O_neg = - sqrt(2/J * (E0-mga*(0.5*phi.^2-1/24*phi.^4)));
phid_5O_pos = + sqrt(2/J * (E0-mga*(0.5*phi.^2-1/24*phi.^4+1/720*phi.^6)));
phid_5O_neg = - sqrt(2/J * (E0-mga*(0.5*phi.^2-1/24*phi.^4+1/720*phi.^6)));
Phasenplan darstellen:
In [6]:
figure('position',[0 0 1000 500])
plot(phi,phid_ex_pos,'k','DisplayName','sin(\varphi)'); hold on;
plot(phi,phid_ex_neg,'k');
plot(phi,phid_1O_pos,'b','DisplayName','\varphi');
plot(phi,phid_1O_neg,'b');
plot(phi,phid_3O_pos,'g','DisplayName','\varphi - 1/6 \varphi^3');
plot(phi,phid_3O_neg,'g');
plot(phi,phid_5O_pos,'r','DisplayName','\varphi - 1/6 \varphi^3 + 1/120 \varphi^5');
plot(phi,phid_5O_neg,'r');
xlabel('\varphi'); ylabel('d \varphi / dt'); title('analytische Loesung')
h = legend("show","location", "eastoutside"); set(h, "interpreter", "tex");
drawnow();
Numerische Lösungen¶
Anfangsbedingungen:
In [7]:
phi0 = 0;
phid0 = sqrt(2*E0/J);
z0 = [phi0; phid0];
ODE Solver Options
In [8]:
lsode_options('absolute tolerance',1e-5);
lsode_options('relative tolerance',1e-5);
lsode_options('integration method','stiff');
lsode_options('initial step size',1e-3);
lsode_options('maximum step size',1e-3);
lsode_options('minimum step size',1e-3);
lsode_options() % aktuelle Konfiguration ausgeben
Zeitvektor
In [9]:
tstart = 0;
tende = 20;
tausgabe = 1e-1;
t = tstart:tausgabe:tende;
Differentialgleichungen als Funktionen $\dot{z} = f(z,t)$.$\quad$ $z=(\varphi,\; \dot{\varphi})^T$,$\quad$ $\dot{z}=(\dot{\varphi},\; \ddot{\varphi})^T$
In [10]:
% exakte DGL
function zd = dgl_ex(z, t, mga, J)
zd = zeros(2,1);
zd(1) = z(2);
zd(2) = -mga/J * sin(z(1));
endfunction
% Naeherung 1. Ordnung
function zd = dgl_1O(z, t, mga, J)
zd = zeros(2,1);
zd(1) = z(2);
zd(2) = -mga/J * z(1);
endfunction
% Naeherung 3. Ordnung
function zd = dgl_3O(z, t, mga, J)
zd = zeros(2,1);
zd(1) = z(2);
zd(2) = -mga/J * (z(1)-1/6*z(1)^3);
endfunction
% Naeherung 5. Ordnung
function zd = dgl_5O(z, t, mga, J)
zd = zeros(2,1);
zd(1) = z(2);
zd(2) = -mga/J * (z(1)-1/6*z(1)^3+1/120*z(1)^5);
endfunction
ODEs lösen:
In [11]:
z_ex = lsode(@(z,t) dgl_ex(z,t,mga,J), z0, t);
z_1O = lsode(@(z,t) dgl_1O(z,t,mga,J), z0, t);
z_3O = lsode(@(z,t) dgl_3O(z,t,mga,J), z0, t);
z_5O = lsode(@(z,t) dgl_5O(z,t,mga,J), z0, t);
Ergebnisse aus der Lösung für $\varphi(t)$ berechnen:
$E_{pot}$ berechnen:
In [12]:
Epot_ex = E0 - 0.5*J*z_ex(:,2).^2;
Epot_1O = E0 - 0.5*J*z_1O(:,2).^2;
Epot_3O = E0 - 0.5*J*z_3O(:,2).^2;
Epot_5O = E0 - 0.5*J*z_5O(:,2).^2;
Lösungen darstellen:
Lösung im Zeitbereich:
In [13]:
figure('position',[0 0 1000 420])
plot(t,z_ex(:,1),'k','DisplayName','sin(\varphi)');
hold on;
plot(t,z_1O(:,1),'b','DisplayName','\varphi');
plot(t,z_3O(:,1),'g','DisplayName','\varphi - 1/6 \varphi^3');
plot(t,z_5O(:,1),'r--','DisplayName','\varphi - 1/6 \varphi^3 + 1/120 \varphi^5');
xlabel('t [s]'); ylabel('\varphi'); title('numerische Loesung');
h = legend("show","location", "eastoutside"); set(h, "interpreter", "tex");
$E_{pot}$
In [14]:
figure('position',[0 0 1000 420])
plot(z_ex(:,1),Epot_ex,'k','DisplayName','sin(\varphi)');
hold on;
plot(z_1O(:,1),Epot_1O,'b','DisplayName','\varphi');
plot(z_3O(:,1),Epot_3O,'g','DisplayName','\varphi - 1/6 \varphi^3');
plot(z_5O(:,1),Epot_5O,'r--','DisplayName','\varphi - 1/6\varphi^3 + 1/120\varphi^5');
xlabel('\varphi'); ylabel('Epot'); title('numerische Loesung');
h = legend("show","location", "eastoutside"); set(h, "interpreter", "tex");
Phasenplan
In [15]:
figure('position',[0 0 1000 420])
plot(z_ex(:,1),z_ex(:,2),'k','DisplayName','sin(\varphi)');
hold on;
plot(z_1O(:,1),z_1O(:,2),'b','DisplayName','\varphi');
plot(z_3O(:,1),z_3O(:,2),'g','DisplayName','\varphi - 1/6 \varphi^3');
plot(z_5O(:,1),z_5O(:,2),'r--','DisplayName','\varphi - 1/6\varphi^3 + 1/120\varphi^5');
xlabel('\varphi'); ylabel('d \varphi / dt'); title('numerische Loesung');
h = legend("show","location", "eastoutside"); set(h, "interpreter", "tex");
