当前位置：首页 > news >正文

Python实例题：Python计算偏微分方程

news 来源：原创 2025/6/15 7:59:21

Python实例题

题目

Python计算偏微分方程

代码实现

import numpy as np
import sympy as sp
from sympy import symbols, Function, diff, Array
from sympy.tensor import tensorcontraction, tensorproduct
from sympy.tensor.array import derive_by_array
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3Dclass TensorAnalysis:"""张量分析计算类，支持张量基本运算、微积分和爱因斯坦求和约定"""def __init__(self):"""初始化张量分析计算器"""pass# ================ 张量基本运算 ================def tensor_product(self, A, B):"""计算两个张量的张量积参数:A: 第一个张量B: 第二个张量返回:numpy数组或sympy数组: 张量积结果"""return tensorproduct(A, B)def tensor_contraction(self, T, indices):"""计算张量的缩并参数:T: 张量indices: 要缩并的指标对列表，每个指标对是一个二元组返回:numpy数组或sympy数组: 缩并结果"""result = Tfor i, j in indices:result = tensorcontraction(result, (i, j))return resultdef tensor_contract_product(self, A, B, contraction_indices):"""计算张量积并进行缩并（张量乘法）参数:A: 第一个张量B: 第二个张量contraction_indices: 要缩并的指标对列表返回:numpy数组或sympy数组: 张量乘法结果"""product = self.tensor_product(A, B)return self.tensor_contraction(product, contraction_indices)def raise_index(self, T, metric):"""指标上升（使用度规张量）参数:T: 张量metric: 度规张量返回:numpy数组或sympy数组: 指标上升后的张量"""# 假设我们提升第一个下指标return self.tensor_contract_product(metric, T, [(1, 2)])def lower_index(self, T, metric_inverse):"""指标下降（使用逆度规张量）参数:T: 张量metric_inverse: 逆度规张量返回:numpy数组或sympy数组: 指标下降后的张量"""# 假设我们降低第一个上指标return self.tensor_contract_product(metric_inverse, T, [(1, 2)])# ================ 张量微积分 ================def tensor_gradient(self, T, coords):"""计算张量的梯度参数:T: 张量coords: 坐标变量列表返回:numpy数组或sympy数组: 张量梯度"""return derive_by_array(T, coords)def tensor_divergence(self, T, coords):"""计算张量的散度参数:T: 张量coords: 坐标变量列表返回:numpy数组或sympy数组: 张量散度"""gradient = self.tensor_gradient(T, coords)# 缩并最后两个指标得到散度return tensorcontraction(gradient, (len(T.shape)-1, len(T.shape)))def tensor_curl(self, vector_field, coords):"""计算向量场的旋度（仅适用于三维空间）参数:vector_field: 向量场coords: 坐标变量列表 (x, y, z)返回:numpy数组或sympy数组: 向量场的旋度"""if len(vector_field) != 3 or len(coords) != 3:raise ValueError("旋度计算仅适用于三维向量场")P, Q, R = vector_fieldx, y, z = coordscurl_x = diff(R, y) - diff(Q, z)curl_y = diff(P, z) - diff(R, x)curl_z = diff(Q, x) - diff(P, y)return Array([curl_x, curl_y, curl_z])def laplacian(self, scalar_field, coords):"""计算标量场的拉普拉斯算子参数:scalar_field: 标量场coords: 坐标变量列表返回:sympy表达式: 标量场的拉普拉斯"""result = 0for coord in coords:result += diff(scalar_field, coord, coord)return result# ================ 爱因斯坦求和约定 ================def einstein_sum(self, expr, *operands):"""使用爱因斯坦求和约定计算张量运算参数:expr: 求和表达式，如"ij,jk->ik"表示矩阵乘法operands: 参与运算的张量列表返回:numpy数组: 计算结果"""return np.einsum(expr, *operands)# ================ 克里斯托费尔符号 ================def christoffel_symbols(self, metric, coords):"""计算克里斯托费尔符号 Γ^k_ij参数:metric: 度规张量 g_ijcoords: 坐标变量列表返回:sympy数组: 克里斯托费尔符号"""n = len(coords)metric_inv = sp.Matrix(metric).inv()# 创建克里斯托费尔符号数组christoffel = Array(np.zeros((n, n, n), dtype=object))# 计算每个分量for k in range(n):for i in range(n):for j in range(n):term = 0for m in range(n):term += 0.5 * metric_inv[k, m] * (diff(metric[i, m], coords[j]) + diff(metric[j, m], coords[i]) - diff(metric[i, j], coords[m]))christoffel[k, i, j] = term.simplify()return christoffel# ================ 黎曼曲率张量 ================def riemann_tensor(self, metric, coords):"""计算黎曼曲率张量 R^l_ijk参数:metric: 度规张量 g_ijcoords: 坐标变量列表返回:sympy数组: 黎曼曲率张量"""n = len(coords)christoffel = self.christoffel_symbols(metric, coords)# 创建黎曼曲率张量数组riemann = Array(np.zeros((n, n, n, n), dtype=object))# 计算每个分量for l in range(n):for i in range(n):for j in range(n):for k in range(n):# 第一项: ∂Γ^l_ik/∂x^jterm1 = diff(christoffel[l, i, k], coords[j])# 第二项: ∂Γ^l_ij/∂x^kterm2 = diff(christoffel[l, i, j], coords[k])# 第三项: Γ^m_ik Γ^l_mjterm3 = 0for m in range(n):term3 += christoffel[m, i, k] * christoffel[l, m, j]# 第四项: Γ^m_ij Γ^l_mkterm4 = 0for m in range(n):term4 += christoffel[m, i, j] * christoffel[l, m, k]# 组合所有项riemann[l, i, j, k] = (term1 - term2 + term3 - term4).simplify()return riemann# ================ 里奇张量和标量曲率 ================def ricci_tensor(self, metric, coords):"""计算里奇张量 R_ij参数:metric: 度规张量 g_ijcoords: 坐标变量列表返回:sympy数组: 里奇张量"""riemann = self.riemann_tensor(metric, coords)# 缩并第一和第三指标得到里奇张量return tensorcontraction(riemann, (0, 2))def scalar_curvature(self, metric, coords):"""计算标量曲率 R参数:metric: 度规张量 g_ijcoords: 坐标变量列表返回:sympy表达式: 标量曲率"""ricci = self.ricci_tensor(metric, coords)metric_inv = sp.Matrix(metric).inv()# 计算 R = g^ij R_ijscalar = 0for i in range(len(coords)):for j in range(len(coords)):scalar += metric_inv[i, j] * ricci[i, j]return scalar.simplify()# ================ 可视化 ================def visualize_vector_field(self, vector_field, coords, domain, density=15):"""可视化二维或三维向量场参数:vector_field: 向量场函数或表达式coords: 坐标变量列表domain: 定义域，格式为[(x_min, x_max), (y_min, y_max), ...]density: 箭头密度"""if len(coords) == 2:# 二维向量场x_min, x_max = domain[0]y_min, y_max = domain[1]x = np.linspace(x_min, x_max, density)y = np.linspace(y_min, y_max, density)X, Y = np.meshgrid(x, y)# 计算向量场值if callable(vector_field):U, V = vector_field(X, Y)else:# 假设vector_field是sympy表达式u_expr, v_expr = vector_fieldu_func = sp.lambdify(coords, u_expr, 'numpy')v_func = sp.lambdify(coords, v_expr, 'numpy')U = u_func(X, Y)V = v_func(X, Y)# 可视化plt.figure(figsize=(10, 8))plt.quiver(X, Y, U, V, np.sqrt(U**2 + V**2), cmap='viridis')plt.colorbar(label='向量大小')plt.xlabel(f'{coords[0]}')plt.ylabel(f'{coords[1]}')plt.title('向量场可视化')plt.grid(True)plt.show()elif len(coords) == 3:# 三维向量场x_min, x_max = domain[0]y_min, y_max = domain[1]z_min, z_max = domain[2]x = np.linspace(x_min, x_max, density)y = np.linspace(y_min, y_max, density)z = np.linspace(z_min, z_max, density)X, Y, Z = np.meshgrid(x, y, z)# 计算向量场值if callable(vector_field):U, V, W = vector_field(X, Y, Z)else:# 假设vector_field是sympy表达式u_expr, v_expr, w_expr = vector_fieldu_func = sp.lambdify(coords, u_expr, 'numpy')v_func = sp.lambdify(coords, v_expr, 'numpy')w_func = sp.lambdify(coords, w_expr, 'numpy')U = u_func(X, Y, Z)V = v_func(X, Y, Z)W = w_func(X, Y, Z)# 可视化fig = plt.figure(figsize=(10, 8))ax = fig.add_subplot(111, projection='3d')ax.quiver(X, Y, Z, U, V, W, length=0.1, normalize=True, cmap='viridis')ax.set_xlabel(f'{coords[0]}')ax.set_ylabel(f'{coords[1]}')ax.set_zlabel(f'{coords[2]}')ax.set_title('三维向量场可视化')plt.show()else:raise ValueError("只能可视化二维或三维向量场")# 示例使用
def example_usage():ta = TensorAnalysis()print("\n===== 张量基本运算示例 =====")# 创建张量A = Array([[1, 2], [3, 4]])B = Array([[5, 6], [7, 8]])print("张量 A:")print(A)print("\n张量 B:")print(B)# 张量积product = ta.tensor_product(A, B)print("\n张量积 A⊗B:")print(product)# 张量缩并contraction = ta.tensor_contraction(product, [(1, 2)])print("\n缩并结果 (A⊗B)_{ik} = A_{ij}B_{jk}:")print(contraction)# 爱因斯坦求和约定einsum_result = ta.einstein_sum('ij,jk->ik', np.array(A), np.array(B))print("\n使用爱因斯坦求和约定 (A⊗B)_{ik}:")print(einsum_result)print("\n===== 张量微积分示例 =====")# 定义坐标变量x, y, z = symbols('x y z')coords = [x, y, z]# 定义标量场scalar_field = x**2 + y**2 + z**2print("\n标量场 φ =", scalar_field)# 计算梯度gradient = ta.tensor_gradient(scalar_field, coords)print("\n梯度 ∇φ =", gradient)# 定义向量场vector_field = Array([x*y, y*z, x*z])print("\n向量场 F =", vector_field)# 计算散度divergence = ta.tensor_divergence(vector_field, coords)print("\n散度 ∇·F =", divergence)# 计算旋度curl = ta.tensor_curl(vector_field, coords)print("\n旋度 ∇×F =", curl)# 计算拉普拉斯算子laplacian = ta.laplacian(scalar_field, coords)print("\n拉普拉斯算子 ∇²φ =", laplacian)print("\n===== 微分几何示例 =====")# 定义二维球面的度规张量theta, phi = symbols('theta phi')metric = Array([[1, 0], [0, sp.sin(theta)**2]])print("\n球面度规张量 g_ij:")print(metric)# 计算克里斯托费尔符号christoffel = ta.christoffel_symbols(metric, [theta, phi])print("\n克里斯托费尔符号 Γ^k_ij:")print(christoffel)# 计算里奇张量ricci = ta.ricci_tensor(metric, [theta, phi])print("\n里奇张量 R_ij:")print(ricci)# 计算标量曲率scalar = ta.scalar_curvature(metric, [theta, phi])print("\n标量曲率 R =", scalar)print("\n===== 向量场可视化示例 =====")# 可视化二维向量场vector_field_2d = Array([-y, x])  # 旋转向量场ta.visualize_vector_field(vector_field_2d, [x, y], [(-3, 3), (-3, 3)], density=10)# 可视化三维向量场vector_field_3d = Array([-x, -y, z])  # 向心+轴向向量场ta.visualize_vector_field(vector_field_3d, [x, y, z], [(-2, 2), (-2, 2), (0, 4)], density=8)if __name__ == "__main__":example_usage()

1. 有限差分法求解偏微分方程

有限差分法是求解 PDE 最基本的方法，它将导数近似为差分，从而将 PDE 转化为代数方程组。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3Ddef solve_heat_equation(L, T, nx, nt, k, initial_condition, boundary_conditions):"""求解一维热传导方程 u_t = k * u_xx参数:L: 空间域长度T: 时间域长度nx: 空间网格点数nt: 时间网格点数k: 热扩散系数initial_condition: 初始条件函数boundary_conditions: 边界条件函数(左边界和右边界)返回:tuple: (x坐标, t坐标, 温度分布u)"""# 网格划分dx = L / (nx - 1)dt = T / nt# 检查稳定性条件if k * dt / dx**2 > 0.5:print(f"警告: 不满足稳定性条件，可能导致数值不稳定")print(f"当前参数: k*dt/dx² = {k*dt/dx**2}，应小于等于0.5")# 创建网格x = np.linspace(0, L, nx)t = np.linspace(0, T, nt+1)# 初始化解数组u = np.zeros((nt+1, nx))# 设置初始条件u[0, :] = initial_condition(x)# 设置边界条件left_bc, right_bc = boundary_conditionsfor j in range(nt+1):u[j, 0] = left_bc(t[j])u[j, -1] = right_bc(t[j])# 显式差分方法求解for j in range(nt):for i in range(1, nx-1):u[j+1, i] = u[j, i] + k * dt / dx**2 * (u[j, i+1] - 2*u[j, i] + u[j, i-1])return x, t, u# 示例：求解热传导方程
def example_heat_equation():# 参数设置L = 1.0         # 空间域长度T = 0.1         # 时间域长度nx = 51         # 空间网格点数nt = 1000       # 时间网格点数k = 0.01        # 热扩散系数# 初始条件：u(x,0) = sin(πx)def initial_condition(x):return np.sin(np.pi * x)# 边界条件：u(0,t) = 0, u(1,t) = 0def left_bc(t):return 0def right_bc(t):return 0boundary_conditions = (left_bc, right_bc)# 求解PDEx, t, u = solve_heat_equation(L, T, nx, nt, k, initial_condition, boundary_conditions)# 可视化结果plt.figure(figsize=(12, 8))# 3D表面图X, T = np.meshgrid(x, t)fig = plt.figure(figsize=(12, 8))ax = fig.add_subplot(111, projection='3d')surf = ax.plot_surface(X, T, u, cmap='viridis')ax.set_xlabel('x')ax.set_ylabel('t')ax.set_zlabel('u(x,t)')ax.set_title('热传导方程的解 u_t = k*u_xx')fig.colorbar(surf, shrink=0.5, aspect=5)plt.show()# 不同时间的温度分布曲线plt.figure(figsize=(10, 6))times = [0, 0.02, 0.05, 0.1]for t_val in times:idx = int(t_val * nt / T)plt.plot(x, u[idx, :], label=f't = {t_val:.2f}')plt.xlabel('x')plt.ylabel('u(x,t)')plt.title('不同时间的温度分布')plt.legend()plt.grid(True)plt.show()if __name__ == "__main__":example_heat_equation()

2. 使用有限元法求解 PDE

对于复杂几何形状和边界条件，有限元法 (FEM) 是更有效的方法。Python 的fenics库提供了强大的有限元求解功能。

from fenics import *
import matplotlib.pyplot as pltdef solve_poisson_equation():"""求解泊松方程 -∇²u = f"""# 创建网格和定义函数空间mesh = UnitSquareMesh(8, 8)  # 8x8的网格V = FunctionSpace(mesh, 'P', 1)  # 一阶拉格朗日有限元# 定义边界条件def boundary(x, on_boundary):return on_boundaryu_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)bc = DirichletBC(V, u_D, boundary)# 定义变分问题u = TrialFunction(V)v = TestFunction(V)f = Constant(-6.0)a = dot(grad(u), grad(v))*dxL = f*v*dx# 计算解u = Function(V)solve(a == L, u, bc)# 计算误差u_e = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)error_L2 = errornorm(u_e, u, 'L2')print(f"L2误差: {error_L2:.6f}")# 计算最大误差vertex_values_u_e = u_e.compute_vertex_values(mesh)vertex_values_u = u.compute_vertex_values(mesh)error_max = np.max(np.abs(vertex_values_u_e - vertex_values_u))print(f"最大误差: {error_max:.6f}")# 可视化结果plt.figure(figsize=(10, 8))plot(u, title='泊松方程的解')plt.colorbar(plot(u))plt.show()# 可视化网格plt.figure(figsize=(8, 6))plot(mesh, title='计算网格')plt.show()if __name__ == "__main__":solve_poisson_equation()

3. 使用谱方法求解 PDE

谱方法适用于周期性边界条件的问题，具有很高的精度。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeintdef solve_kdv_equation():"""求解Korteweg-de Vries方程 u_t + u*u_x + u_xxx = 0"""# 空间域和网格L = 20N = 256dx = L / Nx = np.arange(0, L, dx)# 初始条件: 孤子波c = 1.0  # 波速u0 = 3*c / (np.cosh(0.5 * np.sqrt(c) * (x - L/2))**2)# 波数k = np.fft.fftfreq(N, dx) * 2*np.pi# 定义ODE系统 (在傅里叶空间)def kdv(u_hat, t, k):# 转换回物理空间u = np.fft.ifft(u_hat)# 计算导数u_x = np.fft.ifft(1j * k * u_hat)# 非线性项nonlin = -u * u_x# 转换回傅里叶空间nonlin_hat = np.fft.fft(nonlin)# 色散项dispersion = 1j * k**3 * u_hat# 演化方程du_hat_dt = nonlin_hat + dispersionreturn du_hat_dt# 初始条件的傅里叶变换u0_hat = np.fft.fft(u0)# 时间点t = np.linspace(0, 10, 50)# 求解ODE系统u_hat_sol = odeint(kdv, u0_hat, t, args=(k,))# 转换回物理空间u_sol = np.zeros((len(t), N))for i in range(len(t)):u_sol[i, :] = np.real(np.fft.ifft(u_hat_sol[i, :]))# 可视化结果plt.figure(figsize=(12, 8))plt.imshow(u_sol, extent=[0, L, 0, 10], aspect='auto', cmap='viridis', origin='lower')plt.colorbar(label='u(x,t)')plt.xlabel('x')plt.ylabel('t')plt.title('KdV方程的解 (孤子演化)')plt.show()# 绘制特定时间的解plt.figure(figsize=(10, 6))times = [0, 2, 4, 6, 8, 10]for i, t_val in enumerate(times):idx = int(t_val * (len(t) - 1) / 10)plt.plot(x, u_sol[idx, :], label=f't = {t_val}')plt.xlabel('x')plt.ylabel('u(x,t)')plt.title('不同时间的KdV方程解')plt.legend()plt.grid(True)plt.show()if __name__ == "__main__":solve_kdv_equation()

4. 使用`scipy.integrate.solve_bvp`求解边值问题

对于二阶 PDE 转换的边值问题，可以使用scipy的 BVP 求解器。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_bvpdef solve_bratu_problem():"""求解Bratu问题 u'' + λ*exp(u) = 0，u(0)=u(1)=0"""# 定义方程def fun(x, y, p):lambda_ = p[0]return np.vstack((y[1], -lambda_ * np.exp(y[0])))# 定义边界条件def bc(ya, yb, p):lambda_ = p[0]return np.array([ya[0], yb[0], ya[1] + lambda_])  # ya[1] + lambda_ 是为了得到非平凡解# 初始网格和猜测解x = np.linspace(0, 1, 10)y_guess = np.zeros((2, x.size))  # y[0]是u，y[1]是u'# 求解BVPlambda_guess = 2.0  # 参数λ的初始猜测值sol = solve_bvp(fun, bc, x, y_guess, p=[lambda_guess], verbose=2)# 打印结果print(f"求解成功: {sol.success}")print(f"参数 λ = {sol.p[0]:.6f}")# 可视化解x_plot = np.linspace(0, 1, 100)y_plot = sol.sol(x_plot)plt.figure(figsize=(10, 6))plt.plot(x_plot, y_plot[0], 'b-', linewidth=2)plt.xlabel('x')plt.ylabel('u(x)')plt.title(f'Bratu问题的解 (λ = {sol.p[0]:.4f})')plt.grid(True)plt.show()if __name__ == "__main__":solve_bratu_problem()

5. 使用神经网络求解 PDE (Physics-Informed Neural Networks)

对于复杂 PDE，还可以使用神经网络求解，即 Physics-Informed Neural Networks (PINNs)。

import numpy as np
import tensorflow as tf
import matplotlib.pyplot as pltclass PINN:"""Physics-Informed Neural Network求解器"""def __init__(self, layers, x_train, t_train, u_train, lb, ub):"""初始化PINN参数:layers: 神经网络层数x_train, t_train, u_train: 训练数据lb, ub: 定义域的下界和上界"""# 网络参数self.layers = layersself.lb = lbself.ub = ub# 训练数据self.x_train = x_trainself.t_train = t_trainself.u_train = u_train# 创建神经网络self.model = self.build_model()# 优化器self.optimizer = tf.keras.optimizers.Adam(learning_rate=0.001)def build_model(self):"""构建神经网络模型"""model = tf.keras.Sequential()# 输入层model.add(tf.keras.layers.InputLayer(input_shape=(2,)))# 归一化层model.add(tf.keras.layers.Lambda(lambda x: 2.0 * (x - self.lb) / (self.ub - self.lb) - 1.0))# 隐藏层for layer_size in self.layers:model.add(tf.keras.layers.Dense(layer_size, activation='tanh',kernel_initializer='glorot_normal'))# 输出层model.add(tf.keras.layers.Dense(1))return modeldef grad(self, x, t):"""计算梯度"""with tf.GradientTape(persistent=True) as tape:tape.watch(x)tape.watch(t)u = self.model(tf.concat([x, t], 1))u_x = tape.gradient(u, x)u_t = tape.gradient(u, t)u_xx = tape.gradient(u_x, x)del tapereturn u, u_t, u_xxdef loss(self, x, t, u):"""计算损失函数"""# 预测解u_pred, u_t, u_xx = self.grad(x, t)# PDE残差f = u_t - 0.01 * u_xx  # 热传导方程 u_t = 0.01*u_xx# 损失函数mse_u = tf.reduce_mean(tf.square(u - u_pred))mse_f = tf.reduce_mean(tf.square(f))return mse_u + mse_fdef train_step(self, x, t, u):"""执行一个训练步骤"""with tf.GradientTape() as tape:loss_value = self.loss(x, t, u)# 计算梯度grads = tape.gradient(loss_value, self.model.trainable_variables)# 应用梯度self.optimizer.apply_gradients(zip(grads, self.model.trainable_variables))return loss_valuedef train(self, epochs):"""训练模型"""for epoch in range(epochs):loss_value = self.train_step(self.x_train, self.t_train, self.u_train)if epoch % 100 == 0:print(f'Epoch {epoch}, Loss: {loss_value.numpy():.8f}')def predict(self, x, t):"""预测解"""u_pred = self.model(tf.concat([x, t], 1))return u_preddef solve_heat_equation_pinn():"""使用PINN求解热传导方程"""# 定义定义域lb = np.array([0.0, 0.0])  # 下界 [x_min, t_min]ub = np.array([1.0, 0.2])  # 上界 [x_max, t_max]# 创建训练数据x_train = np.linspace(lb[0], ub[0], 100)[:, None]t_train = np.linspace(lb[1], ub[1], 50)[:, None]X, T = np.meshgrid(x_train, t_train)x_train = X.flatten()[:, None]t_train = T.flatten()[:, None]# 初始条件: u(x,0) = sin(πx)u_train = np.sin(np.pi * x_train) * np.exp(-np.pi**2 * t_train)# 定义神经网络结构layers = [20, 20, 20, 20, 20]  # 隐藏层# 创建并训练PINNpinn = PINN(layers, x_train, t_train, u_train, lb, ub)pinn.train(epochs=5000)# 预测和解的可视化x_plot = np.linspace(lb[0], ub[0], 200)t_plot = np.linspace(lb[1], ub[1], 100)X_plot, T_plot = np.meshgrid(x_plot, t_plot)x_flat = X_plot.flatten()[:, None]t_flat = T_plot.flatten()[:, None]u_pred = pinn.predict(x_flat, t_flat)U_pred = u_pred.numpy().reshape(X_plot.shape)# 可视化结果plt.figure(figsize=(12, 8))plt.imshow(U_pred, extent=[lb[0], ub[0], lb[1], ub[1]], aspect='auto', cmap='viridis', origin='lower')plt.colorbar(label='u(x,t)')plt.xlabel('x')plt.ylabel('t')plt.title('使用PINN求解热传导方程')plt.show()# 特定时间的解plt.figure(figsize=(10, 6))times = [0.0, 0.05, 0.1, 0.15, 0.2]for t_val in times:t_idx = np.abs(t_plot - t_val).argmin()plt.plot(x_plot, U_pred[t_idx, :], label=f't = {t_val}')plt.xlabel('x')plt.ylabel('u(x,t)')plt.title('不同时间的热传导方程解')plt.legend()plt.grid(True)plt.show()if __name__ == "__main__":solve_heat_equation_pinn()