3 Δυναμικά συστήματα 3.3 Μη γραμμικά συστήματα 4 Δυναμική Φορτίων - Δυναμική Δινών

3.4 Μελέτη

3.4.1 Ασκήσεις

Άσκηση 3.1.

Τυχαίο γραμμικό σύστημα
Θεωρήστε το γραμμικό σύστημα

\dot{x}_{1}=ax_{1}+bx_{2}\,,\qquad\dot{x}_{2}=cx_{1}+dx_{2}.

Κληρώστε 4 ακέραιους αριθμούς από το -10 έως το 10 και θεωρήστε ότι οι αριθμοί που κληρώθηκαν δίνουν τις τιμές των $a, b, c, d$ . Βρείτε το είδος και την ευστάθεια του σημείου ισορροπίας για το γραμμικό σύστημα που κληρώσατε. Σχεδιάστε, με μολύβι, με όσο καλύτερη ακρίβεια μπορείτε το διάγραμμα φάσεων (δείξτε με σαφήνεια τον ασταθή και ευσταθή υπόχωρο).

Άσκηση 3.2.

Διαγωνοποίηση πίνακα γραμμικού συστήματος - σάγμα
Για τα γραμμικό σύστημα του παραδείγματος 3.5 σχεδιάστε το διάγραμμα φάσεων στις μεταβλητές $(x_{1},x_{2})$ και επίσης στις $(y_{1},y_{2})$ .

Άσκηση 3.3.

Πίνακας συστήματος για κίνηση φορτίου σε μαγνητικό πεδίο
Το σύστημα Εξ. (4.6) στο Κεφ. 4 δίνει την ταχύτητα φορτίου σε μαγνητικό πεδίο. Βρείτε τη γενική λύση μελετώντας τον πίνακα του συστήματος. Κάνετε το ίδιο για το σύστημα Εξ. (4.47) όπου περιλαμβάνεται και τριβή.

Άσκηση 3.4.

Γραμμικό σύστημα
Μελετήστε το ακόλουθο γραμμικό σύστημα

\dot{x}_{1}=x_{1}+x_{2},\qquad\dot{x}_{2}=-x_{1}+3x_{2}.

Βρείτε τους αναλοίωτους υποχώρους. Σχεδιάστε το διάγραμμα φάσεων.

Άσκηση 3.5.

Μη-γραμμικό σύστημα με τετραγωνικό όρο
Βρείτε τα σημεία ισορροπίας και σχεδιάστε το διάγραμμα φάσεων του μη γραμμικού συστήματος

\dot{x}_{1}=x_{1}x_{2}-4,\qquad\dot{x}_{2}=x_{1}-x_{2}.

3.4.2 Εργασίες

Εργασία 3.1.

Γραμμικό σύστημα
Για το ακόλουθο γραμμικό σύστημα

\dot{x}=-x-3y,\qquad\dot{y}=2y,

(α) Βρείτε τα σημεία ισορροπίας, (β) βρείτε τους αναλλοίωτους υποχώρους και (γ) σχεδιάστε το διάγραμμα φάσεων λύνοντας αριθμητικά το σύστημα.

Βιβλιογραφία: [15].

Εργασία 3.2.

Αναρμονικός ταλαντωτής
Οι παρακάτω εξίσωση περιγράφει έναν αναρμονικό ταλαντωτή:

\ddot{x}=-x+\epsilon x^{3},\qquad\epsilon>0.

(α) Βρείτε τα σημεία ισορροπίας και μελετήστε την ευστάθειά τους. (β) Γύρω από ένα από τα ασταθή σημεία ισορροπίας σχεδιάστε τον ευσταθή και ασταθή υπόχωρο του γραμμικού συστήματος. (γ) Σχεδιάστε το διάγραμμα φάσεων. (δ) Περιγράψτε (και ποσοτικά) τη σχέση του συστήματος με τον αρμονικό ταλαντωτή και με το απλό εκκρεμές.

[Υπόδειξη: Κάνετε όλα τα σχέδια με υπολογιστή. Χρησιμοποιήστε μία τιμή για το $\epsilon$ της επιλογής σας.]

Βιβλιογραφία: [17].

3.4.3 Αλγόριθμοι / Κώδικες

% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
% MATLAB code
% - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

% - - - - - - - - - - - - - - - - - - - -
% solution of a linear system of two equations
% plot the phase portrait
% - - - - - - - - - - - - - - - - - - - -

% define the set of time instances for plotting data points
times = 0:0.1:5;

hold on

% use a set of initial conditions
for x0=-3:3
[time, xy] = ode45('node_func',times,[x0 3]);
plot(xy(:,1),xy(:,2),'-b')
end

for x0=-4:3
[time, xy] = ode45('node_func',times,[x0 -3]);
plot(xy(:,1),xy(:,2),'-b')
end

% format the graph
xmin = -3; xmax = 3;
ymin = -3; ymax = 3;
axis([xmin xmax ymin ymax]);
axis equal;
box on;
xlabel('x','FontSize',18);
ylabel('y','FontSize',18);
set(gca,'XTick',xmin:1:xmax,'FontSize',18);
set(gca,'YTick',ymin:1:ymax,'FontSize',18);

hold off

% print the graph
%print('-depsc','dynsys')


% function for a linear system of equations
function rhs = node_func( t, xx )

x1 = xx(1);
x2 = xx(2);

xdot = x1-2*x2;
ydot = 3*x1-4*x2;

rhs = [xdot; ydot];

end


# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# python code
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
"""
find solution for a dynamical system (2x2)
and render phase portrait
Plot lines in real-time: ask for initial condition and plot - repeatedly
"""

from scipy.integrate import odeint
from numpy import array,linspace
import matplotlib.pylab as plt                # for plotting commands


# system of two ODEs
def deriv(y,t): # return derivatives of the array y
    return array([ 2+y[0]*(1.0-y[1]), y[0]-y[1] ])


# Main program

# prepare graph
plt.figure(figsize=(6,6))
plt.xlabel(r'$x$',fontsize=20,labelpad=0)
plt.ylabel(r'$y$',fontsize=20,labelpad=0)
xmin = -4.0; xmax = 4.0
ymin = -4.0; ymax = 4.0
dx = 2.0; dy = 2.0
plt.xticks(plt.arange(xmin,xmax+dx,dx),fontsize=14)
plt.yticks(plt.arange(ymin,ymax+dy,dy),fontsize=14)

plt.axis([xmin,xmax,ymin,ymax])
plt.ion()
plt.show()

while True:

# input initial condition
    init = input('Give initial condition x,y (0,0 to end): ')
    if init[0] == 0.0 and init[1] == 0.0: break
    yinit = array(init)             # initial values
    #time_interval = input('Give time interval t0,t1 (same to end): ')
    #if time_interval[0] == time_interval[1]: break
    time_interval = (0.0,5.0)
    ntime = 100    # number of points on time interval
    time = linspace(time_interval[0],time_interval[1],ntime)

# integrate ODEs
    y = odeint(deriv,yinit,time)
    print "Final point: ", y[ntime-1,0],",",y[ntime-1,1]

# discard points ourside plotting interval
    yp1 = []
    yp2 = []
    for iline in range(ntime):
        if abs(y[iline,0]) > 1.e4 or abs(y[iline,1]) > 1.e4:
            break
        else:
            yp1.append(y[iline,0])
            yp2.append(y[iline,1])

# plot line
    plt.plot(yp1,yp2,'-b')
    plt.pause(0.01)

#savefig('plot1.png', dpi=96)