# Schémas numériques
# ------------------
# Il est plus simple de définir d'abord une version scalaire des schémas (c'est à dire un schéma pour approcher une
# EDO scalaire) :
def euler_exp_s(f,t0,tf,y0,h):
T=np.append(np.arange(t0,tf,h),tf)
N=T.size
y=np.zeros(N)
y[0]=y0
for n in range(N-1):
y[n+1]=y[n]+dt*f(y[n],T[n])
return T,y
# méthode d'EE programmée vectoriellement, pour un système d'EDO. Elle peut être aussi utilisée dans le cas scalaire
# si on définit y0 comme un numpy array ( par exemple, si y0=1 il faut alors définir y0=np.array([1.]) ).
def euler_exp(f,t0,tf,y0,h):
T=np.append(np.arange(t0,tf,h),tf)
N=T.size
y=np.zeros((N,y0.size))
y[0,:]=y0
for n in range(N-1):
y[n+1,:]=y[n,:]+h*f(y[n,:],T[n])
return T,y
# REMARQUE : Si on a un système de 2 équations, par exemple, de la forme | x'=f1(t,x,y),
# | y'=f2(t,x,y),
# si on appelle X=(x,y), la solution approchée est (X_0,...,X_N)= ( (x_0,...,x_N) , (y_0,...,y_N) )
# Dans la fonction euler_exp qu'on a construit, cette solution approchée est gardée dans le tableau y qui
# a N lignes et 2 colonnes (comme dans odeint).
# Donc dans la colonne 1 il y a (x_0,...,x_N), les valeurs approchées de la première composante de X,
# dans la colonne 2 il y a (y_0,...,y_NN), les valeurs approchées de la deuxième composante de X.
def Heun(f,t0,tf,y0,h):
T=np.append(np.arange(t0,tf,h),tf)
N=T.size
y=np.zeros((N,y0.size))
y[0,:]=y0
for n in range(N-1):
p1=f(y[n,:],T[n])
p2=f(y[n,:]+h*p1,T[n+1])
y[n+1,:]=y[n,:]+h*(p1+p2)/2
return T,y
def RK4(f,t0,tf,y0,h):
T=np.append(np.arange(t0,tf,h),tf)
N=T.size
y=np.zeros((N,y0.size))
y[0,:]=y0
for n in range(N-1):
p1=f(y[n,:],T[n])
p2=f(y[n,:]+h*p1/2,T[n]+h/2)
p3=f(y[n,:]+h*p2/2,T[n]+h/2)
p4=f(y[n,:]+h*p3,T[n]+h)
y[n+1,:]=y[n,:]+(p1+2*p2+2*p3+p4)/6
y[n+1]=y[n]+h*(p1+2*p2+2*p3+p4)/6
return T,y
# Fonctions f second membre de l'EDO
# ----------------------------------
# Q2 - (P1)
def fS1(y,t,b=2.,c=1.):
return c*y*(1.-y/b)
# Q3 - (P2)
def fS3(y,t): # Ecriture vectorielle u=y, v=y'
return np.array([y[1],-y[0]+np.cos(t)])
# Solutions exactes des EDO
# -------------------------
# Q2
def u1(t,a=0.1,b=2.,c=1.):
return b/(1+(np.exp(-c*t))*(b-a)/a)
# Q3
def u3(t):
return np.sin(t)*t/2+5*np.cos(t)+np.sin(t)
# Paramètres Ex 5 - Q2
# --------------------
t0=0
tf=15
a=0.1
h=0.2
y0=np.array([a])
# Solutions approchées par EE, Heun et RK4
# ----------------------------------------
T1E,y1E=euler_exp(fS1,t0,tf,y0,h)
T2E,y2E=euler_exp(fS1,t0,tf,y0,h/2)
T1H,y1H=Heun(fS1,t0,tf,y0,h)
T2H,y2H=Heun(fS1,t0,tf,y0,h/2)
T1RK4,y1RK4=RK4(fS1,t0,tf,y0,h)
T2RK4,y2RK4=RK4(fS1,t0,tf,y0,h/2)
# pour représenter la solution exacte
pas=10**(-4)
vect=np.arange(t0,tf+pas,pas)
plt.figure(1)
plt.plot(vect,u1(vect),label='sol ex')
plt.plot(T1E, y1E,'g',label='h='+str(h))
plt.plot(T2E, y2E,'r',label='h/2='+str(h/2))
plt.legend()
plt.title('euler explicite pour equation logistique')
plt.figure(2)
plt.plot(vect,u1(vect),label='sol ex')
plt.plot(T1H, y1H,'g',label='h='+str(h))
plt.plot(T2H, y2H,'r',label='h/2='+str(h/2))
plt.legend()
plt.title('Heun pour equation logistique')
plt.figure(3)
plt.plot(vect,u1(vect),label='sol ex')
plt.plot(T1RK4, y1RK4,'g',label='h='+str(h))
plt.plot(T2RK4, y2RK4,'r',label='h/2='+str(h/2))
plt.legend()
plt.title('RK4 pour equation logistique')
