Info Mp* : Révisions Maths II Centrale

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as intgr
 
def euler(a, b, N) :
    h = (b - a) / N
    x = np.zeros(N + 1)
    y = np.zeros(N + 1)
    yp = np.ones(N + 1)
    x[0] = a
    yp[0] = 1
    for k in range(N) :
        x[k+1]  = a + (k+1)*h
        yp[k+1] = yp[k] + y[k]*h/(1-x[k])**3
        y[k+1]  = y[k] + yp[k]*h
    plt.plot(x,y)
    # directement avec odeint
    yo = intgr.odeint(lambda Y, x : [Y[1], Y[0]/(1-x)**3], [a, b], x)
    plt.plot(x,yo[:,0]) # valeurs de y
    plt.show()

import numpy as np
import numpy.linalg as alg
 
def J(n) :
    return np.matrix([[1 if j == (n - 1) - i else 0 for j in range(n)] for i in range(n)])
 
def randMatrix(n, p) :
    return np.matrix([[np.random.randint(100) for i in range(p)] for i in range(n)])
 
def centro(A) :
    n,p = A.shape
    assert n == p, "La matrice doit être carrée"
    Jn = J(n)
    return Jn * A * Jn
 
def est_sym(a) :
    return np.allclose(a, a.transpose())
 
def est_anti_sym(a) :
    return np.allclose(a, - a.transpose())
 
def est_cent_sym(a) :
    return np.allclose(a, centro(a))
 
def est_cent_anti(a) :
    return np.allclose(a, -centro(a))
 
 
def decompo(A) :
    tA = A.transpose()
    sA = (A + tA) / 2
    aA = (A - tA) / 2
    d = [(sA + centro(sA)) / 2, (sA - centro(sA))/2, (aA + centro(aA)) / 2, (aA - centro(aA))/2]
    assert est_sym(d[0]) and est_cent_sym(d[0]), "première composante :-("
    assert est_sym(d[1]) and est_cent_anti(d[1]), "seconde composante :-("
    assert est_anti_sym(d[2]) and est_cent_sym(d[2]), "troisième composante :-("
    assert est_anti_sym(d[3]) and est_cent_anti(d[3]), "quatrième composante :-("
    assert np.allclose(sum(d), A), "somme :-("
    return d
 
def blocs(A,B,C,D) :
    n = A.shape[0]
    def m(i,j) :
        if i < n and j < n :
            return A[i, j]
        if i < n and j >= n :
            return B[i, j - n]
        if i >= n  and j < n :
            return C[i - n, j]
        return D[i - n, j - n]
    return np.matrix([[m(i,j) for j in range(2*n)] for i in range(2*n)])
 
def Q(n) :
    In = np.matrix([[1 if i == j else 0 for j in range(n)] for i in range(n)])
    Jn = J(n)
    return blocs(In,-Jn,Jn,In)

import numpy.random as rd
 
def f1a(p, k) :
    S = 0
    while S < k :
        S += 1 + rd.binomial(1, p)
    return S == k
 
def test(n, p, k) :
    print("1 / E(Y1) = " + str(1 / (1 + p)))
    return sum([f1a(p, k) for _ in range(n)]) / n

import numpy.random as rd
import numpy as np
 
def partie1(n, N, S) :
    lancer = rd.binomial(1, 0.5, N)
    prevs = [rd.binomial(1, 0.5, N) for joueur in range(n)]
    reussis = [sum(1 - abs(prev - lancer)) for prev in prevs]
    best = max(reussis)
    gagnants = [score == best for score in reussis]
    gain = S / sum(gagnants)
    return np.array([gagnant * gain for gagnant in gagnants])
 
def simul1(n, N, S, parties) :
    return [sum([partie1(n, N, S) for partie in range(parties)]) / parties]