######################################
#   Thème 2 : Méthodes numériques    #
######################################

## Importations
import numpy as np
import math as m
import matplotlib.pyplot as plt

## I) Méthode des rectangles
# exercices 23/24

# avec le nombre de rectangles en argument :
def rectangle(f,a,b,n) :
    '''donne une approximation de l'integrale de f sur [a,b] par la methode des rectangles en utilisant n rectangles.'''
    S = 0
    h = (b-a)/n
    for k in range(n) :
        S += f(a+k*h)
    return h*S

def f1(x) :
    return 2-x**2

def f2(x) :
    return np.log(x)-1

def f3(x) :
    return x**0.5 - 2

def f4(x) :
    return x-2*np.sin(2*x)

def f5(x) :
    return 1/(1+x**2)

n = 100
print('Avec ',n,' rectangles :')
print("I1 = ",rectangle(f1,0,2,n))
print("I2 = ",rectangle(f2,1,2,n))
print("I3 = ",rectangle(f3,1,6,n))
print("I4 = ",rectangle(f4,0,3,n))
print("I5 = ",rectangle(f5,0,1,n))
print('\n')

# avec le pas en argument :
def rectangle2(f,a,b,pas) :
    '''donne une approximation de l'integrale de f sur [a,b] par la methode des rectangles en utilisant le pas en argument.'''
    S = 0
    n = int((b-a)/pas)
    for k in range(n) :
        S += f(a+k*pas)
    return pas*S

h = 10**-3
print('Avec un pas de ',h,' :')
print("I1 = ",rectangle2(f1,0,2,h))
print("I2 = ",rectangle2(f2,1,2,h))
print("I3 = ",rectangle2(f3,1,6,h))
print("I4 = ",rectangle2(f4,0,3,h))
print("I5 = ",rectangle2(f5,0,1,h))

# exercices 25
def f(n) :
    def fn(x) :
        return x**2*np.sin(n*x)
    return fn

def u(n) :
    return rectangle(f(n),0,1,100)

def liste(n) :
    L=[]
    for k in range(n+1) :
        L.append(u(k))
    return L

# on peut conjecturer que lim u_n = 0
# Ce résultat se démontre par une IPP
# en dérivant x -> x² et en intégrant
# x -> sin(nx).
# (voir : théorème de Riemann-Lebesgue)

## II) Fonctions définies par une intégrale

# exercice 26

def F1(x) :
    def g1(t) :
        return np.exp(-t**2)
    return rectangle(g1,0,x,100)

def F2(x) :
    def g2(t) :
        return 1/(1+t**2+t**4)
    return rectangle(g2,0,x,100)

def trace(f,a,b,N) :
    '''trace la courbe de f sur l'intervalle [a,b] en utilisant N points. '''
    X = np.linspace(a,b,N)
    Y = [f(x) for x in X]
    plt.plot(X,Y)
    plt.show()

# trace(F1,0,1,100)
# trace(F2,-5,5,100)

# exercice 27

def F3(x) :
    def g3(t) :
        return np.sin(4*t)
    return rectangle(g3,x,x**2,1000)

def sol_exacte(x) :
    return (np.cos(4*x)-np.cos(4*x**2))/4

trace(F3,0,10,1000)
trace(sol_exacte,0,10,1000)

# il faut augmenter le nombre de rectangles
# pour que la solution approchée soit précise,
# car la solution exacte a des variations rapides.

## exercice 28

def H1(x) :
    def h1(t) :
        return np.exp(-t*x)/(1+t)
    return rectangle(h1,0,1,100)

def H2(x) :
    def h2(t) :
        return np.sin(x+t)/(x+t)
    return rectangle(h2,0,1,100)

def H3(x) :
    def h3(t) :
        return np.exp(-x*(1+t**2))/(1+t**2)
    return rectangle(h3,0,1,100)

# trace(H1,1,5,1000)
# trace(H2,2,10,1000)
trace(H3,-1,1,1000)

## III) Méthode d'Euler

# exercice 30
def Euler(ED,a,b,y0,pas) :
    '''Resout l'equation differentielle 'ED' par la methode d'Euler sur [a,b], avec y(a) = y0 et en utilisant le pas. On renvoie un couple d'une liste des abscisses t et une liste des ordonnees y(t).'''
    T = [a]
    Y = [y0]
    n = int((b-a)/pas)
    for _ in range(n) :
        t,y = T[-1],Y[-1]
        yp = ED(t,y)
        T.append(t+pas)
        Y.append(y+yp*pas)
    return T,Y

def EDa(t,y) :
    return y

def EDb(t,y) :
    return 2*y

def EDc(t,y) :
    return -y+t

def EDd(t,y) :
    return -t*y+t*np.exp(-t**2/2)/2

def EDe(t,y) :
    '''utilise une fonction z, variable globale.'''
    return -2*y+z(t)

def courbe_integrale(ED,a,b,y0,h) :
    T,Y = Euler(ED,a,b,y0,h)
    plt.plot(T,Y)
    plt.show()

# exercice 31

def systeme_differentiel(SD,a,b,x0,y0,pas) :
    T = [a]
    X = [x0]
    Y = [y0]
    n = int((b-a)/pas)
    for _ in range(n) :
        t = T[-1]
        x = X[-1]
        y = Y[-1]
        xp,yp = SD(x,y)
        T.append(t+pas)
        X.append(x+pas*xp)
        Y.append(y+pas*yp)
    return T,X,Y

def SD(x,y) :
    xp = -2*x+y
    yp = 2*x-2*y
    return xp,yp

def courbe_parametrique(SD,a,b,x0,y0,pas) :
    T,X,Y = systeme_differentiel(SD,a,b,x0,y0,pas)
    plt.plot(X,Y)
    plt.show()

# exercice 32

h = 10**-2
# courbe_integrale(EDa,0,1,1,h)
# courbe_integrale(EDb,0,1,1,h)
# courbe_integrale(EDc,0,1,1,h)
# courbe_integrale(EDd,0,1,1,h)
# courbe_parametrique(SD,0,1,1,0,h)

## exercice 33

def ED33(t,y) :
    return -2*t*y+t*m.exp(-t**2)*y**3

def f33(x) :
    return (3/(m.exp(-x**2)+2*m.exp(2*x**2)))**0.5

T,Y = Euler(ED33,0,2,1,2*10**-3)
Z = [f33(t) for t in T]
plt.plot(T,Y,label = 'solution approchée')
plt.plot(T,Z,'--',label = 'solution exacte')
plt.legend()
plt.show()

## exercice 34

def ED34(t,y) :
    return -2*y/t+t**2

def f34(x) :
    return (x**-2+x**3)/5

h = 5
b = 200
T,Y = Euler(ED34,1,b,2/5,h)
Z = [f34(t) for t in T]
plt.plot(T,Y,label = 'solution approchée')
plt.plot(T,Z,'--',label = 'solution exacte')
plt.legend()
plt.show()

# on pourrait utiliser un pas négatif
# dans la méthode d'Euler pour construire
# une solution approchée en partant d'une
# condition initiale y(c) = y0 avec a < c < b.

## exercice 35

def SD35(x,y) :
    xp = -2*x+y/4
    yp = -3*x
    return xp,yp

a,b,h = 0,10,0.1
def x35(t) :
    return -3*m.exp(-t/2)/4 + 7*m.exp(-3*t/2)/4
def y35(t) :
    return -9*m.exp(-t/2)/2 + 7*m.exp(-3*t/2)/2
T,X,Y = systeme_differentiel(SD35,a,b,1,-1,h)
plt.plot(T,X,label = 'courbe approchée de x')
plt.plot(T,Y,label = 'courbe approchée de y')
Xex = [x35(t) for t in T]
Yex = [y35(t) for t in T]
plt.plot(T,Xex,'--',label = 'courbe exacte de x')
plt.plot(T,Yex,'--',label = 'courbe exacte de y')
plt.legend()
plt.show()