# -*- coding: utf-8 -*-
"""

@author: le prince
"""
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integr


"""Premier exemple"""

T = 30.0 ; t0 = 0 ; y0 = 1.0 ; n = 500
h = T/n   

# Une première méthode sans numpy et sans F1
t = t0
liste_t = [t]
y = y0
liste_y = [y0]
for i in range(n):
    y = y + h * (np.sin(t)-y)
    t = t+h
    liste_t.append(t)
    liste_y.append(y)
    
plt.plot(liste_t , liste_y)
    



# Avec  numpy et fonction F1 associée au système : 

def F1(y,t):
    return -y+np.sin(t)


temps = np.linspace(t0,t0+T,n+1)
Sol = np.zeros(n+1)
Sol[0] = y0
for i in range(n):
    Sol[i+1] = Sol[i] + h*F1( Sol[i],temps[i])
plt.plot(temps,Sol,'b')

#solution exacte
SolEx = [1.5*np.exp(-t)+0.5*(np.sin(t)-np.cos(t)) for t in temps]
plt.plot(temps,SolEx,'r')

# avec odeint
SolScipy = integr.odeint(F1,y0,temps)
plt.plot(temps,SolScipy)



"""Dimension 2"""



        
        
   

def EulerDim2(t0,T,Y0,n,F):
    h=T/n
    temps = np.linspace(t0 , t0+T , n+1)
    list_Y = np.zeros((2,n+1))
    list_Y[:,0] = Y0
    for i in range(n):
        list_Y[:,i+1] = list_Y[:,i] + h*F(list_Y[:,i], temps[i])
    return list_Y
 


    
def F2(Y,t):
    """avec multiplication matricielle.
    Il faut que Y soit un tableau numpy
    (à la bonne dimension)
    """
    A = np.array([[1,1],[0,-1]])
    return np.dot(A,Y)
"""
#remarques :
A = np.array([[1,1],[0,-1]])
print(np.shape(A))
Z = np.array([[0],[1]])
print(np.shape(Z))
print(np.dot(A,Z))
Y0=np.array([0,1])
print(np.shape(Y0))
print(np.dot(A,Y0))
Y1 = [0,1]
A1=[[1,1],[0,-1]]
prod = np.dot(A1,Y1)
print(prod)
print(type(prod))"""

    
##test
T = 3.
t0 = 0
Y0 = np.array([0,1])
n = 10

Sol = EulerDim2(t0,T,Y0,n,F2)
#print(np.shape(Sol))
#print(np.shape(Sol[:,0]))

temps = np.linspace(0 , t0+T , n+1)
X = Sol[0,:] 
Y = Sol[1,:]  

plt.plot(X,Y,'r')

# Une   méthode sans numpy et sans fonction associée au système

T=3.
t0=0
Y0=[0,1]
n=10
h=T/n

t = t0
temps = [t]
y1 , y2 = Y0
liste_y1 = [y1]
liste_y2 = [y2]
for i in range(n):
    y1 , y2 = y1 + h * (y1+y2) , y2 + h * (-y2)
    t = t+h
    temps.append(t)
    liste_y1.append(y1)
    liste_y2.append(y2)
    
plt.plot(liste_y1 , liste_y2)

    

"""Second ordre

résolution de l'équation y''+4y=0"""


def F3(Y,t):
    return np.array([Y[1],-4*Y[0]])


T=10.
t0=0
Y0=np.array([0.5,0.1])
n=100

Sol2 = EulerDim2(t0,T,Y0,n,F3)

temps = np.linspace(0 , t0+T , n+1)
listY = Sol2[0,:] 
listYprime = Sol2[1,:]  

plt.plot(temps,listY,'r')

plt.plot(listY,listYprime,'r')