Le sujet SIP-TD12-13.pdf

Éléments de corrigé

import random
import math
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import matplotlib.colors as colors

def approx_pi(n) :
    inside = 0
    total = 0
    while total < n :
        x = random.uniform(-1, 1)
        y = random.uniform(-1, 1)
        total += 1
        if math.hypot(x, y) < 1 :
            inside += 1
    almost_pi = 4 * (inside / total)
    rel_error = math.fabs(almost_pi - math.pi) / math.pi
    return (almost_pi, rel_error)

# approx_pi(1000)
# approx_pi(100000)

def plot_monte_carlo() :
    fig = plt.figure()         # Create a figure
    fig.patch.set_facecolor('white') # Set its background to white
    axes = plt.axes()          # Create the axes of the figure
    axes.set_aspect('equal')   # Set their aspect ratio to 1
    axes.set_axis_off()        # Do not display the axes
    axes.set_xlim(left=-1, right=1.1)  # The .1 is for avoiding to clip the rectangle
    axes.set_ylim(bottom=-1.1, top=1)
    # Add the rectangle and the circle
    axes.add_patch(patches.Circle((0, 0), 1, linewidth=1, fill=False, edgecolor='blue'))
    axes.add_patch(patches.Rectangle((-1, -1), 2, 2, linewidth=1, fill=False))

def approx_pi_plot(n) :
    plot_monte_carlo()
    inside = 0
    total = 0
    while total < n :
        x = random.uniform(-1, 1)
        y = random.uniform(-1, 1)
        total += 1
        if math.hypot(x, y) < 1 :
            inside += 1
            plt.plot(x, y, 'ro', markersize=2)
        else:
            plt.plot(x, y, 'go', markersize=2)
        plt.pause(0.05)  # Pause for animation
    plt.show()
    almost_pi = 4 * (inside / total)
    rel_error = math.fabs(almost_pi - math.pi) / math.pi
    return (almost_pi, rel_error)

##
# Mandelbrot set
##
def module_real_img(z) :
    return math.hypot(z.real, z.imag)

def module_conj(z) :
    return math.sqrt((z.conjugate()*z).real)

def mandelfunc(c, z) :
    return z*z + c

def mandelbrot(c, n_max) :
    z = 0
    n = 0
    while n < n_max and abs(z) < 2 :
        n += 1
        z = mandelfunc(c, z)
    return n

def mandelbrot_set(xmin, xmax, ymin, ymax, xnum, ynum, n_max) :
    x = [xmin + i * (xmax - xmin)/xnum for i in range(xnum+1)]
    y = [ymin + i * (ymax - ymin)/ynum for i in range(ynum+1)]
    z = [[mandelbrot(complex(x[i], y[j]), n_max) for i in range(xnum)] for j in range(ynum)]
    return z

def plot_mandelbrot(xmin, xmax, ymin, ymax, xnum, ynum, n_max) :
    z = mandelbrot_set(xmin, xmax, ymin, ymax, xnum, ynum, n_max)
    mycolors = [(1,0,0), (0,1,0), (0,0,1), (0,0,0)]
    colmap = colors.LinearSegmentedColormap.from_list("mandelbrot", mycolors, 16)
    plt.imshow(z, cmap=colmap, interpolation='bicubic')
    plt.show()

# plot_mandelbrot(-2.25, 0.75, -1.25, 1.25, 1500, 1250, 64)

##
# PageRank
##
def subdiv_one(n):
    elems = [random.random() for _  in range(n)]
    norm = sum(elems)
    return [elem / norm for elem in elems]

import numpy as np
def create_transition_matrix(size) :
    # First, we create a matrix whose lines sum up to 1
    m = np.matrix([subdiv_one(size) for i in range(size)])
    # and we transpose it to have columns that sum up to 1
    return m.transpose()

# https://en.wikipedia.org/wiki/Matrix_norm
import numpy.linalg as la

def page_rank(size, epsilon) :
    m = create_transition_matrix(size)
    s0 = [[1/size] for _ in range(size)]  # 1 column matrix
    rank = 0
    norm = la.norm((m ** (rank + 1) - m ** rank), np.inf)
    while (norm/size > epsilon) :
        rank += 1
        norm = la.norm((m ** (rank + 1) - m ** rank), np.inf)
    return (m ** rank) * s0

# More efficient version, which does not recompute m**k at
# each iteration. It also compare the norm directly to epsilon,
# without dividing by the number of pages
def page_rank_lowcost(size, epsilon) :
    m = create_transition_matrix(size)
    mk = m
    s0 = [[1/size] for _ in range(size)]
    rank = 0
    norm = la.norm(m, np.inf)
    while (norm > epsilon) :
        rank += 1
        mk2 = mk * m
        norm = la.norm((mk2 - mk), np.inf)
        mk = mk2
    return mk * s0

##
# Prey-predator system
##
# The function that gives the derivative from the current values
def f(x, y, a, b) :
    return x * (a - b * y)

# The function that gives the next value of z (which is either x or y)
# from the current value of z, of the other variable (t), the parameters
# of the problem (a and b) and the integration step (Euler explicit scheme).
def next_step(zn, tn, a, b, step) :
    return zn + step * f(zn, tn, a, b)

def populations(x0, y0, alpha, beta, delta, gamma, step, n) :
    x, y = x0, y0
    lx, ly = [x0], [y0]
    for _ in range(n) :
        # We have dx/dt =  alpha*x - beta*x*y
        # and     dy/dt = -gamma*y + delta*y*x
        # hence the parameters of next_step
        x, y = next_step(x, y, alpha, beta, step), \
               next_step(y, x, -gamma, -delta, step)
        lx.append(x)
        ly.append(y)
    return lx, ly

def display_populations() :
    # Initial populations
    x0, y0 = 0.9, 0.9
    # Parameters
    alpha, beta, delta, gamma = 2/3, 4/3, 1, 1
    # Temporal resolution
    resolution = 1000
    # Compute until that time
    time_max = 20
    # Time step size
    step = 1 / resolution
    # Total number of steps
    n = time_max * resolution
    ## Compute the populations
    prey, pred = populations(x0, y0, alpha, beta, delta, gamma, step, n)
    # Plot the result
    time = [step * i for i in range(n+1)]
    plt.subplot(121)
    # First subplot for the evolution of populations through time
    plt.title("Population versus time")
    plt.plot(time, prey, label="Preys")
    plt.plot(time, pred, label="Predators")
    plt.xlabel("Time")
    plt.ylabel("Population")
    plt.legend(loc='upper left')
    # Second subplot for the phase diagram
    plt.subplot(122)
    plt.title("Predators versus preys")
    plt.plot(prey, pred)
    plt.xlabel("Preys")
    plt.ylabel("Predators")
    plt.show()