# ======================================================================= #
# ======================================================================= #
# Résolution numérique d'un système différentiel linéaire dont la matrice #
# associée est une matrice réelle carrée d'ordre 2.                       #
# Version de base.                                                        #
# 06/12/2014                                                                 #
# ======================================================================= #
# ======================================================================= #

# Importation des bibliothèques requises.
# =======================================
import numpy as np
import matplotlib.pyplot as plt

# La fonction donnant la dérivée x' en fonction de x et de y.
# ===========================================================
def fx(x,y):
    global a, b
    return a * x + b * y

# La fonction donnant la dérivée y' en fonction de x et de y.
# ===========================================================
def fy(x,y):
    global c, d
    return c * x + d * y

# Le coeur de l'algorithme : le schéma numérique proprement dit.
# ==============================================================
def Euler2D(x0,y0,n,pas):
    # Initialisation des tableaux des abscisses et des ordonnées (n+1 valeurs car il y a n pas ...).
    # Ces tableaux sont initialisés avec des 0.
    x = np.zeros(n+1)
    y = np.zeros(n+1)
    x[0], y[0] = x0, y0
    # Le schéma numérique.
    for i in range(n):
        x[i+1] = x[i] + pas * fx(x[i],y[i])
        y[i+1] = y[i] + pas * fy(x[i],y[i])
    plt.plot(x,y,'white')

# ==================== #
# Programme principal. #
# ==================== #

# Initialisations
# ===============
# Les coefficients de la matrice A et les conditions initiales
global a, b, c, d
# Les valeurs propres sont réelles strictement positives.
# Coefficients de la matrice A
#a, b = 7, -2
#c, d = 15, -4
# Conditions initiales
#x_ci = [1,1,-1]
#y_ci = [0,3,1]

# Les valeurs propres sont réelles non nulles de signes contraires.
# Coefficients de la matrice A
#a, b = 2, 1
#c, d = -3, -2
# Conditions initiales
#x_ci = [1,1,-2]
#y_ci = [0,-1,1]

# Les valeurs propres sont réelles, l'une d'elle est nulle.
# Coefficients de la matrice A
#a, b = 4, -2
#c, d = 2, -1
# Conditions initiales
#x_ci = [1,1,1]
#y_ci = [0,2,-1]

# Les valeurs propres sont complexes conjuguées de partie réelle non nulle.
# Trajectoire = spirales.
# Coefficients de la matrice A
a, b = 3**0.5/2, -0.5
c, d = 0.5, 3**0.5/2
# Conditions initiales
x_ci = [1,1,1]
y_ci = [0,2,-1]

# Les valeurs propres sont des imaginaires purs conjugués.
# Trajectoires = ellipses.
# Coefficients de la matrice A
#a, b = 0, -5
#c, d = 5, 0
# Conditions initiales
#x_ci = [1,2,5]
#y_ci = [0,0,0]

# Le nombre de points de chaque branche.
n = int(input('Veuillez saisir le nombre de points souhaités : '))
# Le pas (typiquement 0.01). On peut (doit ? :) fournir un pas négatif.
pas = 0
while pas == 0:
    pas = float(input('Veuillez saisir le pas : '))

# On efface le graphique précédent ...
plt.clf()

# Premiers éléments graphiques.
plt.axes(axisbg='black')
plt.grid(c='grey', lw='2', ls='-')

# Génération des branches des trois courbes intégrales.
# =====================================================
for i in range(3):
    Euler2D(x_ci[i],y_ci[i],n,pas)

# AFFICHAGE.
# ==========

# Génération des éléments graphiques généraux (titre, légendes, ...).
plt.xlabel('x')
plt.ylabel('y')
plt.title('Méthode d\'Euler avec '+str(n)+' pas',backgroundcolor = 'w',style = 'italic')
plt.legend(loc = 'lower left',ncol = 2)

# Affichage des éléments graphiques.
plt.show()