🖥️ Oscillateur de Wien ★★

Pour rappel, un oscillateur de Wien est constitué d’un montage amplificateur non-inverseur et d’un filtre de Wien.

On note

• $𝑢\left(𝑡\right)$ la tension en entrée de l’amplificateur non-inverseur (et donc en sortie du filtre de Wien)
• $𝑣\left(𝑡\right)$ la tension en sortie de l’amplificateur non-inverseur (et donc en entrée du filtre de Wien)

Les fonctions de transfert du filtre de Wien et de l’amplificateur non-inverseur sont respectivement données par :

On pose $𝑤=\frac{\mathrm{d}𝑢}{\mathrm{d}𝑡}$. On cherche à mettre le problème sous la forme d’un problème d’Euler de la forme $\frac{\mathrm{d}\overrightarrow{𝑌}}{\mathrm{d}𝑡}=𝐹\left(𝑡,\overrightarrow{𝑌}\right)$ avec $\overrightarrow{𝑌}=\left(\begin{array}{c}𝑢\\ 𝑣\\ 𝑤\end{array}\right)$.

1/ Exprimer les dérivées de $𝑢$, $𝑣$ et $𝑤$ uniquement en fonction de $𝑢$, $𝑣$, $𝑤$, $𝑅$, $𝐶$, ${𝑅}_1$, ${𝑅}_2$, ${𝐴}_0$ et $𝜏$.

Avec Python, on représente le vecteur $\overrightarrow{𝑌}$ par un array numpy à 3 éléments.

2/ Écrire une fonction Python qui prend en entrée Y et renvoie sa dérivée. On pourra supposer les variables $𝑅$, $𝐶$, ${𝑅}_1$, ${𝑅}_2$, ${𝐴}_0$ et $𝜏$ déjà définies.

import numpy as np
def F(t, Y):
  """
  t : le temps (inutile ici, mais solve_ivp le fournit)
  Y : array numpy à 3 éléments contenant u, v et w
  Renvoie un array numpy contenant la dérivée de Y par rapport au temps.
  """
  ...

import numpy as np
def F(t, Y):
  u,v,w = Y
  return np.array([
        w,
        A0/tau * u - 1/((1+R2/R1)*tau/A0) * v,
        (-1/(R*C)**2+1/(R*C)*A0/tau) * u - 1/(R*C)/((1+R2/R1)*tau/A0) * v - 3/(R*C) * w
    ])

Dans un premier temps, on utilise la fonction scipy.integrate.solve_ivp qui résout numériquement des équations différentielles mises sous la forme d’un problème d’Euler. Cette fonction utilise des variantes de la méthode d’Euler la rendant plus précise.

La fonction solve_ivp prend en argument

• la fonction F définie précédemment
• le couple $\left({𝑡}_{𝑖},{𝑡}_{𝑓}\right)$ correspondant à l’intervalle de temps sur lequel on souhaite résoudre l’équation différentielle
• le tableau $\overrightarrow{𝑌}\left({𝑡}_{𝑖}\right)$ correspondant aux conditions initiales

La fonction solve_ivp renvoie un objet. Si on le stocke dans la variable solution,

• solution.t contient les temps
• solution.y contient les valeurs successives de $\overrightarrow{𝑌}$ à ces temps

3/ Définir et affecter les variables $𝑅=1\hspace{0.17em}\text{kΩ}$, $𝐶=10\hspace{0.17em}\text{nF}$, ${𝑅}_1$, ${𝑅}_2$, ${𝐴}_0=\text{100\hspace{0.17em}000}$ et $𝜏=0\text{,01}\hspace{0.17em}\text{s}$ avec des valeurs vraisemblables satisfaisant la condition de démarrage des oscillations.

R = ... # Ohm
C = ... # F
R1 = ... # Ohm
R2 = ... # Ohm
A0 = ...
tau = ... # s

R = 1000 # Ohm
C = 10E-9 # F
R1 = 1000 # Ohm
R2 = 2100 # Ohm, soit un gain de 3,1 : au-dessus du seuil avec de la marge
A0 = 1E5
tau = 0.01 # s
Vsat = 15 # V, tension de saturation de l'ALI

4/ Définir Y0 avec de très petites valeurs pour $𝑢$, $𝑣$ et $𝑤$ ($0\text{,0001}$ par exemple).

Y0 = ...

Y0 = np.array([0.0001, 0.0001, 0.0001])

5/ Définir tf pour observer une dizaine d’oscillations.

tf = ... # s

tf = 10 * 2 * np.pi * R * C

6/ Tracer $𝑢$ et $𝑣$ en fonction du temps.

import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

solution = solve_ivp(..., ..., ...)
plt.plot(..., ..., label='u(t)')
plt.plot(..., ..., label='v(t)')
plt.xlabel('Temps (s)')
plt.ylabel('Tension (V)')
plt.legend()
plt.show()

import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

solution = solve_ivp(F, (0,tf), Y0)
plt.plot(solution.t, solution.y[0], label='u(t)')
plt.plot(solution.t, solution.y[1], label='v(t)')
plt.xlabel('Temps (s)')
plt.ylabel('Tension (V)')
plt.legend()
plt.show()

7/ Vérifier la condition de démarrage des oscillations.

8/ Vérifier la valeur de la période des oscillations.

Dans la suite, on souhaite se passer de la fonction solve_ivp et implémenter nous-même la méthode d’Euler.

On note ${𝑌}_{𝑖}=𝑌\left(𝑖\cdot \mathrm{\Delta }𝑡\right)$ où $\mathrm{\Delta }𝑡$ est la durée entre deux échantillons (période d’échantillonnage).

9/ Dans le cas général, exprimer ${𝑌}_{𝑖+1}$ en fonction de ${𝑌}_{𝑖}$, de $𝑖$, $\mathrm{\Delta }𝑡$ et de la fonction $𝐹$.

10/ Implémenter la méthode d’Euler pour simuler l’évolution des tensions pour un oscillateur de Wien. $\mathrm{\Delta }𝑡$ sera choisi de sorte qu’il y ait environ $200$ échantillons par période.

Delta_t = ... # s, environ 200 échantillons par période
N = int(tf / Delta_t) # Nombre d'échantillons
t = np.zeros(N)
Y = np.zeros((N,3))
t[0] = 0
Y[0] = Y0
for i in range(1,N):
    Y[i] = ...
    t[i] = ...

Delta_t = (2 * np.pi * R * C) / 200
N = int(tf / Delta_t)
t = np.zeros(N)
Y = np.zeros((N,3))
t[0] = 0
Y[0] = Y0
for i in range(1,N):
    Y[i] = Y[i-1] + Delta_t * F(t[i-1], Y[i-1])
    t[i] = t[i-1] + Delta_t

11/ Adapter le code précédent pour prendre en compte la saturation de l’ALI.

Vsat = 15 # V, tension de saturation de l'ALI

def F_sat(t, Y):
  u,v,w = Y
  dv = A0/tau * u - 1/((1+R2/R1)*tau/A0) * v
  if (v >= Vsat and dv > 0) or (v <= -Vsat and dv < 0):
      dv = 0 # la sortie de l'ALI est bloquée, elle ne varie plus
  return np.array([w, dv, -u/(R*C)**2 - 3/(R*C) * w + dv/(R*C)])

tf = 60 * 2 * np.pi * R * C
Delta_t = (2 * np.pi * R * C) / 200
N = int(tf / Delta_t)
t = np.zeros(N)
Y = np.zeros((N,3))
t[0] = 0
Y[0] = Y0
for i in range(1,N):
    Y[i] = Y[i-1] + Delta_t * F_sat(t[i-1], Y[i-1])
    Y[i][1] = min(max(Y[i][1], -Vsat), Vsat)
    t[i] = t[i-1] + Delta_t