######################################################################
# ZÉROS - DICHOTOMIE
######################################################################

def f(x):
    return x**3/2 - 3*x**2 + x + 4

(a, b, N) = (0, 4, 10)

for _ in range(N):
    m = (a+b)/2

    if(f(a)*f(m)<0):
        b = m
    else:
        a = m

print((a+b)/2, "est une valeur approchée d'un zéro à", (b-a)/2**(N+1), "près.")


######################################################################
# ZÉROS - MÉTHODE DE NEWTON-RAPHSON : RACINE CUBIQUE D'UN RÉEL POSITIF
######################################################################

c = 3250

x = 0

while(x**3 < c):
    x += 1

for k in range(5):
    x = x - (x**3-c) / (3*x**2)

print(f"La racine cubique de {c} est environ {x}.")


######################################################################
# CALCUL DIFFÉRENTIEL - MÉTHODE D'EULER EXPLICITE : Y' = -Y + COS(T), Y(0) = 0
######################################################################

from math import *
import matplotlib.pyplot as plt
import numpy as np

(a, b, N) = (0, 2*pi, 10)

def f(t, y):             # y' = f(t, y)
    return -y + cos(t)

y = 0                    # condition initiale y0 en a
Y = [y]
T = np.linspace(a, b, N+1)
h = (b-a) / N

for k in range(0,N):     # Schéma d'Euler explicite
    y = y + h * f(T[k], y)
    Y.append(y)

plt.plot(T, Y, 'o-', ms=80/N, label=f"{N} intervalles")
plt.legend()
plt.show()


######################################################################
# CALCUL DIFFÉRENTIEL - MÉTHODE DE RUNGE-KUTTA 4 : Y' = Y, Y(0) = 1, CALCUL DE E^2 (CODE COMPLET)
######################################################################

from math import *
import numpy as np

(a, b, N) = (0, 2, 100)
u = 1                       # condition initiale y0
T = np.linspace(a, b, N+1)
h = (b-a)/N

def f(t, y):                # y' = f(t, y)
    return y

def RK(t, u):               # Schéma de Runge-Kutta 4
    k1 = f(t, u)
    k2 = f(t+h/2, u+h*k1/2)
    k3 = f(t+h/2, u+h*k2/2)
    k4 = f(t+h, u+h*k3)
    return h/6*(k1 + 2*k2 + 2*k3 + k4)

for k in range(0, N):
    u = u + RK(T[k], u)

print("Valeur connue :  ", e**2)
print("Valeur calculée :", u, "avec", N, "intervalles")
print("Erreur :", e**2-u)


######################################################################
# CALCUL DIFFÉRENTIEL - COMPARAISON AVEC EULER EXPLICITE (REMPLACER LA LIGNE 20 DU CODE PRÉCÉDENT)
######################################################################

u = u + h*f(T[k], u)


######################################################################
# CALCUL INTÉGRAL - MÉTHODE DES RECTANGLES : CALCUL
######################################################################

from math import *

def f(x):
    return 3 + x - sqrt(1+x**3)

a, b = 0, 3         # intervalle
N, d = 100, 1/2    # nb de rectangles, décalage

h = (b-a)/N         # pas
S = 0               # initialisation de la somme

for k in range(N):
    S += f(a+(k+d)*h)

print(h*S)


######################################################################
# CALCUL INTÉGRAL - MÉTHODE DES RECTANGLES : VISUALISATION (CODE COMPLET)
######################################################################

import numpy as np
import matplotlib.pyplot as plt
from math import *

def f(x):
    return 3 - sqrt(1+x**3) + x

a, b = 0, 3
N, d = 10, 1/2

h = (b-a)/N

# Configuration du graphique
plt.grid(lw=0.5)
plt.axhline(y=0, c='black')
plt.axvline(x=0, c='black')

# Tracé des rectangles
for k in range(N):
    RX = [[a+k*h], [a+k*h], [a+(k+1)*h], [a+(k+1)*h]]
    RY = [0, f(a+(k+d)*h), f(a+(k+d)*h), 0]
    plt.plot(RX, RY, 'r', lw='1')

# Tracé de la courbe
X = np.linspace(a, b, 101)
Y = [f(x) for x in X]
plt.plot(X, Y, lw='2')

plt.show()


######################################################################
# CALCUL INTÉGRAL - MÉTHODE DE SIMPSON : CALCUL
######################################################################

from math import *

def f(x):
    return 3 - sqrt(1+x**3) + x

a, b, N = 0, 3, 100

h = (b-a)/N
S = f(a) + f(b)

for k in range(1, N):
    S += 2*f(a + k*h)

for k in range(N):
    S += 4*f(a + (k+1/2)*h)

print(S*h/6)


######################################################################
# SOLUTIONS - EULER EXPLICITE : COMPARAISON DE PLUSIEURS PAS AVEC LA SOLUTION EXACTE
######################################################################

from math import *
import matplotlib.pyplot as plt
import numpy as np

(a, b) = (0, 2*pi)     # intervalle

def f(t, y):        # y' = f(t, y)
    return -y + cos(t)

def s(t):       # fonction solution
    return 1/2*(cos(t) + sin(t) - e**(-t))

for N in [10, 50, 100, 250]:
    y = s(a)      # condition initiale y0 en a
    Y = [y]
    T = np.linspace(a, b, N+1)
    h = (b-a) / N   # pas

    for k in range(0,N):    # Schéma d'Euler expl
        y = y + h * f(T[k], y)
        Y.append(y)

    print(N,"intervalles : erreur =", round(100*(s(b)-Y[-1])/s(b),2), "%")
    plt.plot(T, Y, ls="--", label=f"{N} intervalles")

# Tracé de la fonction solution
L = np.linspace(a, b, 200)
plt.plot(L, [s(t) for t in L], c="black", lw="2")

plt.legend()
plt.show()