Python求解一维热传导方程解析解遇问题，寻求调试帮助

作为Python新手，我正在尝试求解一维热传导方程的解析解，编写了对应的Python代码，但运行结果不符合预期，无法自行排查修复。

给定的解析解公式 #

本人编写的Python代码 #

import numpy as np
import matplotlib.pyplot as plt

print("1D heat equation solver for a wall (Analytical Solution)")

wall_thickness = 31  # Thickness of the wall in centimeters
L = wall_thickness  # Length of the wall
max_time = 1800  # Total simulation time in seconds
delta_x = 0.05
delta_t = 0.01
D = 93 #Diffusivity coeficient

# Define the spatial grid
x = np.arange(0, L + delta_x, delta_x)

# Initialize solution: the grid of u(x, t)
u = np.empty((max_time, len(x)))

# Set the boundary conditions
u_sur = 149  # Surrounding temperature
u_in = 38.0  # Initial temperature inside the wall

# Set the initial condition
u[0] = u_in

Tbase = u_sur + 2 * (u_in - u_sur)
def analytical_solution(x, t):
    T = np.zeros_like(x)
    for m in range(1, 100):  # Sum from m=1 to a large number (e.g., 100)
        T += np.exp(-D*((m * np.pi) / L)**2*t) * ((1 - (-1) ** m) / (m * np.pi)) * np.sin((m * np.pi * x) / L)
    return Tbase + T

# Calculate the analytical solution for each time step
for k in range(1, max_time):
    t = k * delta_t
    u[k] = analytical_solution(x, t)
    #print(f"Valeurs de la sol analytique {u[k]}.\n")

# Plot the temperature profiles at specific time points
plt.figure(figsize=(10, 6))
time_points = [0, 360, 720, 1080, 1440, 1800]
for t in time_points:
    k = int(t / delta_t)
    if k < max_time:
        plt.plot(x, u[k], label=f"t = {t} seconds")

plt.title("Temperature Profiles (Analytical Solution) at Specific Time Points")
plt.xlabel("x (thickness)")
plt.ylabel("Temperature")
plt.ylim(0,150)
plt.xlim(0, L)
plt.legend()
plt.show()

print("Done!")

运行结果 #

可以看到结果不符合预期，图中红色曲线为预期效果，恳请各位提供帮助，谢谢！

内容的提问来源于stack exchange，提问作者Yovan Bach