#### TP n°0: Méthodes numériques ####
#################

### OUTILS GRAPHIQUES ###


## Question 1 
import matplotlib.pyplot as plt    
tab_x=[1,1.1,1.2,1.3,1.5,1.6,1.7,1.8,1.9,2.0]
tab_y=[2,2.1,2.2,2.4,2.8,3.1,3.5,3.8,4.0,4.1]
plt.figure()
plt.scatter(tab_x,tab_y,color="red",marker="*")
plt.show()

## Question 2 
import matplotlib.pyplot as plt
import numpy as np
tab_x=np.linspace(-2,2,1000)
tab_y=tab_x**3
plt.figure()
plt.plot(tab_x,tab_y,color="green")
plt.show()
## Question 3 
import matplotlib.pyplot as plt
import numpy as np
tab_t=np.linspace(0,10e-3,5000)
tab_e=3*np.cos(2000*tab_t)
tab_s=2.2*np.cos(2000*tab_t+0.4)
plt.figure()
plt.plot(tab_t,tab_e,label="e(t)")
plt.plot(tab_t,tab_s,label="s(t)")
plt.legend()
plt.show()
## Question 4 
import matplotlib.pyplot as plt
import numpy as np
Rc,wc=3e-3,9e9
tab_t=np.linspace(0,2*np.pi/wc,1000)
tab_x=Rc*np.sin(wc*tab_t)
tab_y=Rc*(np.cos(wc*tab_t)-1)
plt.figure()
plt.plot(tab_x,tab_y)
plt.axis("equal")
plt.show()

### ÉQUATIONS ALGÉBRIQUES ###


## Question 5 
import matplotlib.pyplot as plt
import numpy as np
def f(x):
    return x*np.exp(x)-2*x**3-3
tab_x=np.linspace(0,3.5,1000)
tab_y=f(tab_x)
plt.figure()
plt.plot(tab_x,tab_y)
plt.show()

a,b=0,3.5
eps=1e-6
while b-a>eps:
    m=(a+b)/2
    if f(a)*f(m)>0:
        a=m
    else:
        b=m
print((a+b)/2)
## Question 6 
import scipy
print(scipy.optimize.bisect(f,0,3.5))
## Question 7 
import matplotlib.pyplot as plt
import numpy as np
def g(x):
    return np.cos(x)-x
tab_x=np.linspace(0,np.pi/2,1000)
tab_y=g(tab_x)
plt.figure()
plt.plot(tab_x,tab_y)
plt.show()
import scipy
print(scipy.optimize.bisect(g,0,np.pi/2))
## Question 8 
import matplotlib.pyplot as plt
import numpy as np
w0,Q=3000,2.1
def H(w):
    return 2/np.sqrt((1-w**2/w0**2)**2+(w/Q/w0)**2)
tab_w=np.linspace(0,10000,10000)
tab_H=H(tab_w)
maxH=np.max(tab_H)
plt.figure()
plt.plot(tab_w,tab_H)
plt.axhline(maxH/np.sqrt(2),linestyle="--")
plt.show()
def f(w):
    return H(w)-maxH/np.sqrt(2)
import scipy
print(scipy.optimize.bisect(f,0,w0))
print(scipy.optimize.bisect(f,w0,10000))
### INTÉGRATION, DÉRIVATION ###


## Question 9 
import math
def f(x):
    return math.sqrt(1-x**2)
integ=0
a,b=0,1
N=1000
dx=(b-a)/N
for i in range(N):
    xi=a+(i+0.5)*dx
    integ=integ+f(xi)
integ=integ*dx
print(integ)
print(math.pi/4)
## Question 10 
a,b=0.1358,0.0000364
T=273.15+25
n,R=1,8.314
Vi,Vf=2.5e-3,20e-3
def P(V):
    return n*R*T/(V-n*b)-a*n**2/V**2
W=0
N=1000
dV=(Vf-Vi)/N
for i in range(N):
    V=Vi+(i+0.5)*dV
    W=W-P(V)
W=W*dV
print(W)

### ÉQUATIONS DIFFÉRENTIELLES ###


## Question 11 
import matplotlib.pyplot as plt
import numpy as np
tau=0.5
pas,tmax=tau/1000,5*tau
N=int(tmax/pas)
tab_t=np.zeros(N+1)
tab_f=np.zeros(N+1)
tab_f[0]=1
for n in range(N):
    tab_t[n+1]=tab_t[n]+pas
    tab_f[n+1]=tab_f[n]+pas*(4-tab_f[n])/tau
plt.figure()
plt.plot(tab_t,tab_f)
plt.xlabel("t")
plt.ylabel("f")
plt.show()
## Question 12 
import matplotlib.pyplot as plt
import numpy as np
pas,tmax=1e-3,10
N=int(tmax/pas)
tab_t=np.zeros(N+1)
tab_f=np.zeros(N+1)
for n in range(N):
    tab_t[n+1]=tab_t[n]+pas
    tab_f[n+1]=tab_f[n]+pas*(1/(1+tab_t[n])-3/2*tab_f[n])
plt.figure()
plt.plot(tab_t,tab_f)
plt.xlabel("t")
plt.ylabel("f")
plt.show()

## Question 13 
import matplotlib.pyplot as plt
import numpy as np
omega0,alpha=2.2,1.2
T0=2*np.pi/omega0
pas,tmax=T0/1000,10
N=int(tmax/pas)
tab_t=np.zeros(N+1)
tab_x=np.zeros(N+1)
tab_xp=np.zeros(N+1)
tab_x[0]=3
for n in range(N):
    tab_t[n+1]=tab_t[n]+pas
    tab_x[n+1]=tab_x[n]+pas*tab_xp[n]
    tab_xp[n+1]=tab_xp[n]+pas*(-omega0**2*tab_x[n]-alpha*tab_xp[n]*abs(tab_xp[n]))
plt.figure()
plt.plot(tab_t,tab_x)
plt.xlabel("t")
plt.ylabel("x")
plt.show()
plt.figure()
plt.plot(tab_x,tab_xp)
plt.xlabel("x")
plt.ylabel("xpoint")
plt.show()

## Question 14 
import matplotlib.pyplot as plt
import numpy as np
import scipy
omega0,alpha=2.2,1.2
T0=2*np.pi/omega0
def derivee(inconnues,t):
    x,xp=inconnues[0],inconnues[1]
    xpp=-omega0**2*x-alpha*xp*abs(xp)
    return [xp,xpp]
ci=[3,0]
tab_t=np.linspace(0,10,10000)
solution=scipy.integrate.odeint(derivee,ci,tab_t)
tab_x=solution[:,0]
tab_xp=solution[:,1]
plt.figure()
plt.plot(tab_t,tab_x)
plt.xlabel("t")
plt.ylabel("x")
plt.show()
plt.figure()
plt.plot(tab_x,tab_xp)
plt.xlabel("x")
plt.ylabel("xpoint")
plt.show()
## Question 15 
import matplotlib.pyplot as plt
import numpy as np
import scipy
omega0,alpha=2.2,1.2
T0=2*np.pi/omega0
def derivee(inconnues,t):
    x,xp=inconnues[0],inconnues[1]
    xpp=-omega0**2*x-alpha*xp*abs(xp)
    return [xp,xpp]
tab_t=np.linspace(0,10,10000)
plt.figure()
for ci in [[3,0],[1.5,2],[-3,0],[-1.5,-2]]:
    solution=scipy.integrate.odeint(derivee,ci,tab_t)
    tab_x=solution[:,0]
    tab_xp=solution[:,1]
    plt.plot(tab_x,tab_xp)
plt.xlabel("x")
plt.ylabel("xpoint")
plt.show()
## Question 16 
import matplotlib.pyplot as plt
import numpy as np
import scipy
omega0,epsilon=2*np.pi,1
def derivee(inconnues,t):
    x,xp=inconnues[0],inconnues[1]
    xpp=epsilon*(1-x**2)*xp-omega0**2*x
    return [xp,xpp]
tab_t=np.linspace(0,5,10000)
co=[0,0]
solution=scipy.integrate.odeint(derivee,ci,tab_t)
tab_x=solution[:,0]
tab_xp=solution[:,1]
plt.figure()
plt.plot(tab_x,tab_xp)
plt.xlabel("x")
plt.ylabel("xpoint")
plt.show()
## Question 17 
import matplotlib.pyplot as plt
import numpy as np
import scipy
omega0=2*np.pi
def derivee(inconnues,t):
    th,thp=inconnues[0],inconnues[1]
    thpp=-omega0**2*np.sin(th)
    return [thp,thpp]
tab_t=np.linspace(0,5,10000)
plt.figure()
for th0 in range(10,171,10):
    ci=[th0*np.pi/180,0]
    solution=scipy.integrate.odeint(derivee,ci,tab_t)
    tab_th=solution[:,0]*180/np.pi
    plt.plot(tab_t,tab_th)
plt.xlabel("t")
plt.ylabel("theta")
plt.show()
### PROBABILITÉS − STATISTIQUES ###


## Question 18 
import matplotlib.pyplot as plt
import numpy as np
tab_x=np.random.normal(3.1,0.1,1000)
plt.figure()
plt.hist(tab_x)
plt.show()
## Question 19 
print(np.mean(tab_x))
print(np.std(tab_x))
## Question 20 
import matplotlib.pyplot as plt
import numpy as np
N=1000
tab_a=np.random.normal(4.1,0.1,N)
tab_b=np.random.normal(3.2,0.2,N)
tab_c=tab_a/tab_b
print(np.mean(tab_c))
print(np.std(tab_c))
plt.figure()
plt.hist(tab_c)
plt.show()
## Question 21 
import numpy as np
N=1000
tab_OA=np.random.normal(-12.5,0.1,N)
tab_OAp=np.random.normal(62.4,0.7,N)
tab_fp=1/(1/tab_OAp-1/tab_OA)
print(np.mean(tab_fp))
print(np.std(tab_fp))

## Question 22 
import matplotlib.pyplot as plt
import numpy as np
data_x=[0.5,1.0,1.5,2.5,3.0,3.5,4.0,4.5,5.0]
data_y=[1.6,2.1,2.3,3.1,3.4,3.5,4.2,4.5,4.7]
coefs=np.polyfit(data_x,data_y,1)
a,b=coefs[0],coefs[1]
print(a,b)
liste_modele=[]
for i in range(len(data_x)):
    liste_modele.append(a*data_x[i]+b)
plt.figure()
plt.scatter(data_x,data_y)
plt.plot(data_x,liste_modele)
plt.show()
## Question 23 
import matplotlib.pyplot as plt
import numpy as np
data_x=[0.5,1.0,1.5,2.5,3.0,3.5,4.0,4.5,5.0]
data_ux=[0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05,0.05]
data_y=[1.6,2.1,2.3,3.1,3.4,3.5,4.2,4.5,4.7]
data_uy=[0.06,0.06,0.09,0.09,0.11,0.12,0.15,0.17,0.20]
liste_a,liste_b=[],[]
N,M=1000,1000    
for m in range(M):
    tab_x=np.random.normal(data_x,data_ux)
    tab_y=np.random.normal(data_y,data_uy)
    coefs=np.polyfit(tab_x,tab_y,1)
    liste_a.append(coefs[0])
    liste_b.append(coefs[1])
print(np.mean(liste_a))
print(np.std(liste_a))
print(np.mean(liste_b))
print(np.std(liste_b))