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

手机工信部网站备案查询数据分析网站html模板下载

wzjs 2025/8/31 8:34:19

手机工信部网站备案查询,数据分析网站html模板下载,淘宝运营培训班学费大概多少,网站建设中目录是什么意思目录 Python实例题题目代码实现 1. 有限差分法求解偏微分方程 2. 使用有限元法求解 PDE 3. 使用谱方法求解 PDE 4. 使用scipy.integrate.solve_bvp求解边值问题 5. 使用神经网络求解 PDE (Physics-Informed Neural Networks) 总结 Python实例题题目 Python计算偏…

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()