import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as spla
import matplotlib.pyplot as plt

L = 1.0 
T = 10 
N = 50 
alpha = 0.01
dt = 0.01 
num_steps = int(T / dt)

h = L / N 
x = np.linspace(0, L, N + 1)

u0 = np.sin(np.pi * x) 

u_left = 0.0
u_right = 0.0

M = np.zeros((N + 1, N + 1))
K = np.zeros((N + 1, N + 1))

for i in range(1, N):
    M[i, i] += h / 3
    M[i, i - 1] += h / 6
    M[i, i + 1] += h / 6

    K[i, i] += 2 * alpha / h
    K[i, i - 1] -= alpha / h
    K[i, i + 1] -= alpha / h

# Assemblage
M = sp.csr_matrix(M)
K = sp.csr_matrix(K)

# Crank-Nicolson
A = M + dt / 2 * K
B = M - dt / 2 * K

# Conditions aux limites (Dirichlet)
A[0, :] = 0
A[0, 0] = 1
A[-1, :] = 0
A[-1, -1] = 1

B[0, :] = 0
B[0, 0] = 1
B[-1, :] = 0
B[-1, -1] = 1

# Résolution temporelle
u = u0.copy()
solution = [u.copy()]

for n in range(num_steps): 
    b = B @ u
    b[0] = u_left
    b[-1] = u_right
    u = spla.spsolve(A, b)
    solution.append(u.copy())

time_points = [n * num_steps // 10 for n in range(11)]

for t in time_points:
    plt.plot(x, solution[t], label=f"t = {(t * dt):.0f}s")
plt.title("Solution Eq Chaleur 1D")
plt.xlabel("x")
plt.ylabel("u(x, t)")
plt.legend()
plt.grid(True)
plt.show()

/Library/Frameworks/Python.framework/Versions/3.12/lib/python3.12/site-packages/scipy/sparse/_index.py:151: SparseEfficiencyWarning: Changing the sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient.
  self._set_arrayXarray(i, j, x)