###Q1###
import matplotlib.pyplot as plt

def f(y, t):
    '''Fonction des deux variables y et t'''
    return y

def Euler1(f,a,t0,T,n):
    '''dans cette version, t est la liste des temps t[i] et y est la liste des ordonnées y[i]'''
    dt = T/n
    t = [t0 + i*dt for i in range(n+1)]
# NB : l'alimentation de t par t += dt dans une boucle for entraînerait des erreurs d'arrondis successifs
    y = [a]
    for i in range(n):
        y.append(y[i] + dt * f(y[i], t[i]))
# NB : on peut aussi écrire "for i in range(1,n+1):" suivi de "y.append(y[i-1]+dt*f(y[i-1],t[i-1]))"
    plt.plot(t,y)
    plt.show()

# ou

def Euler1B(f,a,t0,T,n):
    '''dans cette version, "temps" est la liste des temps t et "ordonnees" est la liste des ordonnées y'''
    dt = T/n
    temps = [t0 + i*dt for i in range(n+1)]
    y = a; ordonnees = [y]
# NB : le calcul qui suit nécessite une première valeur pour y
    for t in temps[:-1]:
# NB: il faut exclure la dernière valeur de temps, non utilisée pour le calcul de yn
        y = y + dt * f(y,t); ordonnees.append(y)
    plt.plot(temps,ordonnees)
    plt.show()

###Q2###

###Q3###
import numpy as np

def suivant(A, B, Xi, h):
    return Xi + h * (np.dot(A, Xi) + B)
# * est la multiplication scalaire, + l'addition matricielle et np.dot la multiplication matricielle de numpy

###Q4###
def Euler(A, B, X0, h, N):
    '''on suppose N >= 1'''
    n = np.shape(A)[0]
    X = np.zeros((n, N+1))
    X[:, 0] = X0
    for i in range(N):
        X[:, i+1] = suivant(A, B, X[:,i], h)
    return X

###Q5###
A = np.array([[3, 4], [5, -2]])
B = np.array([3, 4])
X0 = np.array([1, -1])
t0, tN = 0, 1 #bornes du temps
N = 100
h = (tN - t0)/N
t = [t0 + i*h for i in range(N+1)]
X = Euler(A, B, X0, h, N)
X1, X2 = X[0], X[1]
plt.plot(t, X1)
plt.plot(t, X2)
plt.show()

###Q6###

###Q7###
def suivant(Ai, Bi, Xi, h):
    return Xi + h * (np.dot(Ai, Xi) + Bi)

def Euler_second_ordre(y0, z0, h, la, lb, lc):
    N = len(la)
    X = np.zeros((2, N+1))
    X[:, 0] = np.array([y0, z0]) # il faut entrer une matrice ligne dans la colonne
    for i in range(N):
        Ai = np.array([[0, 1], [-lb[i], -la[i]]])
        Bi = np.array([0, lc[i]])
        X[:, i+1] = suivant(Ai, Bi, X[:, i], h)
    return X

###Q7.1###
# Exemple d'équation d'Euler:
# (t**2)*y'' + (3*t)*y' + y = cos(t)- t * sin(t)
# dont la solution avec y0 = 1 + sin(1) et z0 = cos(1) - sin(1)
# est sin(t)/t + 1/t + ln(t)/t sur R+.
# On écrit : y'' +(3/t) * y' + (1/t**2) y = (cos(t)- t * sin(t)))/t**2
from numpy import cos, sin, log
t0, tN = 1, 30
N = 1000
h = (tN - t0)/N
temps = [t0 + i*h for i in range(N+1)]
la = [3/t for t in temps[:-1]]
lb = [1/t**2 for t in temps[:-1]]
lc = [(cos(t)-t*sin(t))/t**2 for t in temps[:-1]]
y0 = 1 + sin(1)
z0 = cos(1) - sin(1)

X = Euler_second_ordre(y0, z0, h, la, lb, lc)
plt.plot(temps, X[0])

Y = [(sin(t))/t+ (1/t) + (log(t))/t for t in temps]
plt.plot(temps, Y)

plt.show()

###Q8###

###Q9###

###Q10###

###Q11###

###Q12###
def initialise_tableau(L, n):
    N = len(L)
    T = np.zeros((n, N), dtype = np.uint16)
    T[0] = L
    return T

###Q13###
def ligne_suivante(T, i):
    N = np.shape(T)[1]

    # relations Rext
    T[i+1, 0] = (4*T[i, 0] + 2*T[i, 1]) / 6
    T[i+1, N-1] = (2*T[i, N-2] + 4*T[i, N-1]) / 6

    # relations Rint
    for j in range(1, N-1):
        T[i+1, j] = (T[i, j-1] + 4*T[i, j] + T[i, j+1] ) / 6
    # procédure : pas de return

###Q14###
def simulation(L, n):
    T = initialise_tableau(L, n)
    for i in range(n-1):
        ligne_suivante(T, i)
    return T

N = 8
L = [0 for _ in range(N)]; L[N//2] = 2**16 - 1
assert (simulation(L, 10) ==
np.array([[    0,     0,     0,     0, 65535,     0,     0,     0],
          [    0,     0,     0, 10922, 43690, 10922,     0,     0],
          [    0,     0,  1820, 14563, 32767, 14563,  1820,     0],
          [    0,   303,  3640, 15473, 26699, 15473,  3640,   606],
          [  101,   808,  5056, 15371, 22957, 15371,  5106,  1617],
          [  336,  1398,  6067, 14916, 20428, 14924,  6235,  2780],
          [  690,  1999,  6763, 14359, 18592, 14393,  7107,  3931],
          [ 1126,  2574,  7235, 13798, 17186, 13878,  7792,  4989],
          [ 1608,  3109,  7552, 13268, 16070, 13415,  8339,  5923],
          [ 2108,  3599,  7764, 12782, 15160, 13011,  8782,  6728]],
      dtype=np.uint16)).all()

###Q15###

###Q15###

###Fin###

###Exercice###
##Q2
def f(X,t): # fonction de l'equa diff X' = f(X,t)
    '''Avec X = [y, y'], on a X' = [y', y''] = [X[1], 1/(1-t)*X[0]]'''
    return [X[1], 1/(1-t)*X[0]]

from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt

N = 100
liste_t = np.linspace(0, 0.95, N)
X0 = np.array([0,1]) # donné dans l'énoncé
sol = odeint(f, X0, liste_t)

plt.figure()
plt.plot(liste_t, sol[:, 0]) # sol est un np.array (100,2) de première colonnes: les valeurs de y, de seconde: les valeurs de y'

liste_t = np.linspace(0, -0.95, N)
X0 = np.array([0,1]) # donné dans l'énoncé
sol = odeint(f, X0, liste_t)
plt.plot(liste_t, sol[:, 0])
plt.show()

##Q3
liste_n = [n for n in range(N+1)]

a, b = 0, 1
liste_a = [a, b]
for n in range(2, N+1):
    a = ((n-2)/n) * b + (1/(n*(n-1))) * a
    liste_a.append(a)
    a, b = b, a

plt.figure()
plt.scatter(liste_n,liste_a, marker = '+') # on constate 0 <= an <= 1
plt.show()

##Q4
def S100(x):
    S = 0
    for i in range(N+1):
        S = S + liste_a[i] * x**i
    return S

liste_S = [S100(x) for x in liste_t]

plt.figure()
plt.plot(liste_t, liste_S)
plt.plot(liste_t, sol[:, 0])
plt.show()

##
# from numpy import cos, sin, log
# t0, tN = 0, -0.95
# N = 100
# h = (tN - t0)/N
# temps = [t0 + i*h for i in range(N+1)]
# la = [0 for t in temps[:-1]]
# lb = [-1/(1-t) for t in temps[:-1]]
# lc = [0 for t in temps[:-1]]
# y0 = 0
# z0 = 1
#
# X = Euler_second_ordre(y0, z0, h, la, lb, lc)
# plt.plot(temps, X[0])
#
# plt.show()