偏微分方程数值解是数学与应用数学、工程力学等专业的重要课程内容, 其教学既要求学生掌握基本理论, 又要求其具备一定的程序实现与结果分析能力.传统教学模式多以MATLAB等国外软件为主要工具^{[1, 2]}, 在数值计算教学中已形成较为成熟的应用范式.但从当前课程实践情况来看, 此类教学方式仍存在两个较为突出的不足: 一方面, 课堂讲授往往侧重公式推导与离散结果说明, 学生对“弱形式如何落到程序实现”这一关键环节把握不够; 另一方面, 实操教学常依赖较多环境配置与案例准备, 在有限课时内较难形成完整、连贯且可复现的训练过程.

随着课程目标由“会推导、会计算”逐步转向“能建模、能实现、能分析”, 偏微分方程数值解教学的关注重点也应由单纯的理论传授进一步转向理论、算法与实验的贯通.尤其在有限元方法等内容的教学中, 学生虽然能够理解弱形式、单元剖分、矩阵组装等基本概念, 但在将离散思想真正转化为可执行程序时, 仍普遍面临“离散过程抽象、编程链条过长、实验反馈滞后”等问题.如何在有限课时内组织起一条由数学模型、离散思想、程序实现到结果分析的完整教学链条, 已成为偏微分方程数值解课程改革中的一个关键问题.

面对国产化替代与本土化教学资源建设的现实需求, 北太天元作为具有自主知识产权的通用科学计算软件, 为教学改革提供了新的实现路径^[3].与传统仅把软件作为数值实验执行工具的做法不同, 若能够将国产科学计算平台进一步嵌入课程设计之中, 则其意义不仅在于完成算例求解, 更在于为课程实操组织提供统一环境和稳定载体.北太天元兼具矩阵运算、程序编写、图形可视化与PDE相关工具支持等功能, 能够较好覆盖偏微分方程数值解课程中的主要实操环节, 从而为构建本土化、低门槛、易推广的教学路径提供条件.

基于此, 本文不再仅将北太天元视为数值实验的实现工具, 而是将其作为组织偏微分方程数值解实操教学的核心载体.本文以一维和二维非稳态热传导方程为案例, 围绕有限元离散、Crank–Nicolson时间推进、程序实现与结果分析, 构建“理论讲授—问题拆解—代码实践—结果验证—反思提升”相衔接的教学流程, 重点讨论国产软件平台如何服务于偏微分方程数值解课程的实操教学改革.

本文的工作不仅关注数值结果本身的正确性, 也关注教学组织方式的改进及其推广价值.具体而言, 一方面, 本文围绕热传导方程构建了兼具典型性与可操作性的教学案例, 使抽象的有限元离散过程能够在具体问题中得到呈现; 另一方面, 依托北太天元平台完成了从网格剖分、矩阵组装到数值可视化的完整教学实现, 在一定程度上降低了课程实操的组织复杂度.在此基础上, 本文进一步总结了一条可迁移、可复制的国产软件融入路径, 以期为数值计算方法、偏微分方程数值解及相关实操课程提供参考, 并推动教学模式由“结果演示”向“过程训练”转变.

热传导方程是描述热量传递过程的核心偏微分方程, 广泛应用于机械工程热设计、材料热处理、建筑热环境模拟等领域^{[4, 5]}.目前, 数值求解热传导方程的常用方法包括有限差分法、有限元法与边界元法等, 其中有限元法因对复杂几何区域的适应性强、离散精度可控等优点, 是数值计算方法课程中偏微分方程模块的重点内容之一^[6].从教学角度看, 热传导方程既具有明确的物理背景, 又能够较好覆盖弱形式建立、空间离散、时间推进和误差分析等关键教学环节, 因而适合作为偏微分方程数值解课程中的典型实操案例.

下面以有限元法求解非稳态热传导方程为例进行演示.有限元法是求解偏微分方程的经典数值解法, 尤其适用于复杂几何区域与高精度需求的场景.其基本思想是将连续的求解区域划分为有限个“单元”, 在每个单元上构造低次多项式, 即使用形函数对真解进行分片插值, 再通过“弱形式”将偏微分方程转化为线性代数方程组, 最终通过组装单元矩阵得到全局系统并求解, 从而获得离散节点上的近似解^[7].在本文中, 这一理论推导过程不仅是数值分析的基础, 也是后续教学实施的核心内容之一.通过将理论框架、离散格式与软件实现对应起来, 可以使学生更清晰地理解“公式如何落地为程序”的完整过程.

下面我们考虑用有限元法求解非稳态热传导方程^[8]

$ \begin{equation} \frac{\partial u(x, t)}{\partial t} = \alpha \Delta u(x, t) + f(x, t), \quad (x, t) \in \Omega \times (0, T], \end{equation} $

(2.1)

其中$ \alpha>0 $为热扩散系数, $ \Delta $为拉普拉斯算子, $ \Omega\subset\mathbb{R}^d\; (d=1, 2) $为有界求解区域, $ f(x, t) $为热源项.同时, 在区域边界$ \Gamma $上满足边界条件$ u(x, t)=g(x, t), \; (x, t)\in\Gamma\times(0, T], $在初始时刻$ t=0 $时, $ x\in\Omega $中满足初始条件$ u(x, 0)=u_0(x). $上述模型既保留了偏微分方程数值解课程中的典型理论要素, 又便于在后续教学中组织离散推导、程序实现与结果验证, 因此适合作为贯穿全文的案例载体.

考虑一维含时变热源的非稳态热传导方程, 其控制方程为方程(2.1).其中热扩散系数$ \alpha = 2.0 $, 时变热源项为$ f(x, t) = cos(t)cos(x) + 2sin(t)cos(x). $待求的温度场$ u(x, t) $定义在空间区间$ \Omega = [0, 1] $上.边界条件为Dirichlet类型, 在$ x=0 $处, $ u(0, t) = sin(t) $; 在$ x=1 $处, $ u(1, t) = sin(t)cos(1) $.初始条件为$ t=0 $时, $ u(x, 0) = 0 $.该问题的解析真解为$ u_{exact}(x, t) = sin(t)cos(x) $.

对于二维情形下的无热源非稳态热传导方程, 控制方程同样为方程(2.1).其中热扩散系数$ \alpha = 1.0 $, 热源项$ f(x, t) \equiv 0 $, $ \nabla^2 = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} $为二维拉普拉斯算子.待求的温度场$ u(x, y, t) $定义在空间区域$ \Omega = [0, 2] \times [0, 2] $上.边界条件为齐次Dirichlet条件, 即在$ \partial\Omega $($ x=0 $、$ x=2 $、$ y=0 $、$ y=2 $)处有$ u(x, y, t) = 0 $.初始条件为$ t=0 $时, $ u(x, y, 0) = sin\left(\frac{\pi x}{2}\right)sin\left(\frac{\pi y}{2}\right) $.该问题的解析真解为$ u_{exact}(x, y, t) = e^{-\frac{\pi^2 \alpha t}{2}} sin\left(\frac{\pi x}{2}\right)sin\left(\frac{\pi y}{2}\right). $

对求解的一维区域$ \Omega=[0, 1] $进行均匀网格剖分, 单元数$ N_{e} $, 节点数$ N_{n}=N_{e}+1 $, 网格步长$ h=\frac{1}{N_{e}} $. 单元$ e $的形函数^[7]为

$ N_{1}(x)=\frac{x_{2}^{e}-x}{h_{e}}, \quad N_{2}(x)=\frac{x-x_{1}^{e}}{h_{e}} , $

其中$ x_{1}^{e}, x_{2}^{e} $为单元$ e $的左右节点坐标, 形函数导数为$ \frac{dN_{1}}{dx}=-\frac{1}{h_{e}}, \frac{dN_{2}}{dx}=\frac{1}{h_{e}} $. 单元的一致质量矩阵和刚度矩阵分别为

$ M_{e}=\frac{h_{e}}{6}\left[\begin{array}{ll} 2 & 1 \\ 1 & 2 \end{array}\right], \quad K_{e}=\frac{1}{h_{e}}\left[\begin{array}{cc} 1 & -1 \\ -1 & 1 \end{array}\right] . $

对于任意时刻$ t $的单元载荷向量采用1点高斯积分(单元中点$ x_{mid}=\frac{x_{1}^{e}+x_{2}^{e}}{2} $)计算可得

$ \begin{equation} F_{e}(t)=\frac{h_{e}}{2}\left[\begin{array}{l} 1 \\ 1 \end{array}\right] f\left(x_{mid}, t\right). \end{equation} $

(2.2)

对求解域$ \Omega=[0, 2]\times[0, 2] $实施矩形网格剖分, 再将每个矩形单元进一步拆分为两个不重叠的三角形单元.设$ x $方向节点数为$ N_x $, $ y $方向节点数为$ N_y $, 则三角形单元总数为$ N_e=2(N_x-1)(N_y-1) $.

针对任意单个三角形单元, 不妨取单元节点编号为$ 1, 2, 3 $, 其对应的全局坐标分别为$ (x_1, y_1) $、$ (x_2, y_2) $和$ (x_3, y_3) $.局部坐标采用三角形自然坐标$ (\xi, \eta) $, 其需满足约束条件$ \xi\geq0 $、$ \eta\geq0 $且$ \xi+\eta\leq1 $.该单元的3个节点在该局部坐标下的坐标分别为节点1$ (\xi, \eta)=(0, 0) $、节点2$ (\xi, \eta)=(1, 0) $、节点3$ (\xi, \eta)=(0, 1) $.

单元上的形函数取一阶拉格朗日形函数, 形式为$ N(\xi, \eta)=\begin{bmatrix}N_1(\xi, \eta)\\N_2(\xi, \eta)\\N_3(\xi, \eta)\end{bmatrix} $, 其中$ N_i(\xi, \eta) $($ i=1, 2, 3 $)为第$ i $个节点对应的形函数, 需满足节点插值条件$ N_i(\text{节点}j)=\delta_{ij} $($ \delta_{ij} $为克罗内克函数, $ j=1, 2, 3 $).结合上述约束条件与插值条件, 可推导得到二维三角形单元的一阶拉格朗日形函数矩阵^[7]

$ \begin{equation} N(\xi, \eta)=\begin{bmatrix}1-\xi-\eta\\\xi\\\eta\end{bmatrix}. \end{equation} $

(2.3)

全局坐标下单元内任意点的全局坐标$ (x, y) $表述如下

$ \begin{cases} x(\xi, \eta)=N_1x_1+N_2x_2+N_3x_3=(1-\xi-\eta)x_1+\xi x_2+\eta x_3, \\ y(\xi, \eta)=N_1y_1+N_2y_2+N_3y_3=(1-\xi-\eta)y_1+\xi y_2+\eta y_3, \end{cases} $

整理得

$ \begin{cases} x=x_1+\xi(x_2-x_1)+\eta(x_3-x_1), \\ y=y_1+\xi(y_2-y_1)+\eta(y_3-y_1). \end{cases} $

雅可比变换对应的雅可比矩阵$ J $如下

$ \begin{equation} J=\frac{\partial(x, y)}{\partial(\xi, \eta)}=\begin{bmatrix} \frac{\partial x}{\partial \xi} & \frac{\partial x}{\partial \eta} \\ \frac{\partial y}{\partial \xi} & \frac{\partial y}{\partial \eta} \end{bmatrix}=\begin{bmatrix} x_2-x_1 & x_3-x_1 \\ y_2-y_1 & y_3-y_1 \end{bmatrix}. \end{equation} $

(2.4)

根据链式法则, 联立式(2.3)和(2.4), 可得形函数对全局坐标$ (x, y) $的偏导数, 即

$ B = \frac{\partial N}{\partial(x, y)} = J^{-1} \frac{\partial N}{\partial(\xi, \eta)}, $

其中$ J $为雅可比矩阵.单元刚度矩阵的形式为

$ K_{e, ij} = \int_{\Omega_e} \nabla N_i \cdot \nabla N_j \, d\Omega. $

由于$ \nabla N = B $表示形函数在全局坐标系下的梯度, 因此$ \nabla N_i \cdot \nabla N_j $对应矩阵$ BB^\text{T} $的元素.在采用固定三角形单元的假设下, $ B $在单元内为常量, 积分结果为$ BB^\text{T} $乘以单元面积$ A_e $, 故单元刚度矩阵可表示为

$ K_e = BB^\text{T} A_e. $

对应的单元质量矩阵采用集中质量格式以简化计算, 取为

$ M_e = \frac{A_e}{3} I_3, $

其中$ A_e $为单元面积, $ I_3 $为三阶单位矩阵.

借助单元-节点的索引映射关系, 我们将单元质量矩阵与单元刚度矩阵逐步叠加, 构建全局质量矩阵$ M $和全局刚度矩阵$ K $.计算过程中采用北太天元的稀疏矩阵结构存储全局矩阵, 可有效降低内存占用.同时, 根据边界施加的齐次Dirichlet条件, 先识别并提取边界节点, 将其值固定为边界条件对应的零值, 进而提取内部节点对应的子矩阵$ M_\mathrm{int} $和$ K_\mathrm{int} $, 仅针对内部节点求解控制方程, 大幅减少计算量.

在完成空间离散之后, 还需进一步建立适合课堂讲授与程序实现的时间推进格式.为兼顾离散精度、稳定性与实现难度, 本文选取Crank–Nicolson格式作为时间离散方法, 以形成较为完整的教学演示链条.我们对时间区间$ [0, T] $进行均匀离散, 得到时间点序列

$ t^n = n \Delta t, \quad n = 0, 1, \dots, N, \quad N \Delta t = T. $

采用Crank-Nicolson格式, 在时间中点$ t = t^n + \Delta t/2 $处对解的一阶时间导数进行中心差分离散, 并对扩散项及热源项在$ t^n $和$ t^{n+1} $时刻取平均近似, 从而构造出稳定的全离散格式.

对于一维含时变热源的非稳态热传导方程, 其离散形式为

$ \begin{equation} \left(M + \frac{\alpha \Delta t}{2} K\right) u^{n+1} = \left(M - \frac{\alpha \Delta t}{2} K\right) u^n + \frac{\Delta t}{2} \left(F^n + F^{n+1}\right), \end{equation} $

(2.5)

其中各符号含义如下

● $ M $: 全局质量矩阵, 由所有单元质量矩阵组装而成;

● $ K $: 全局刚度矩阵, 由所有单元刚度矩阵组装而成;

● $ F^n $: $ t^n $时刻的全局载荷向量, 表示热源项的空间离散;

● $ u^n $, $ u^{n+1} $: 分别为$ t^n $和$ t^{n+1} $时刻的全局节点温度向量, 存储所有空间节点上的温度值;

● $ \alpha $: 热扩散系数;

● $ \Delta t $: 时间步长.

接下来进一步说明一维情形下热源项的离散与全局载荷向量的组装过程.考虑一维线性两节点单元$ e $, 其端点坐标为$ x_1^e $和$ x_2^e $, 单元长度$ h_e = x_2^e - x_1^e $.该单元上的形函数定义为

$ N_1^e(x) = \frac{x_2^e - x}{h_e}, \quad N_2^e(x) = \frac{x - x_1^e}{h_e}. $

则$ t^n $时刻的单元载荷向量为

$ F_e^n = \int_{x_1^e}^{x_2^e} \begin{bmatrix} N_1^e(x) \\ N_2^e(x) \end{bmatrix} f(x, t^n) dx. $

为简化计算, 对任意时刻$ t $采用中点积分近似取单元中点$ x_{\mathrm{mid}} = \frac{x_1^e + x_2^e}{2} $处的热源值, 并假设其在整个单元内恒定, 则有

$ F_e(t) \approx \frac{h_e}{2} f(x_{\mathrm{mid}}, t) \begin{bmatrix} 1 \\ 1 \end{bmatrix}, $

此即等效集中载荷形式, 与公式(2.3)中的处理方式一致.

最终, 通过将每个单元的局部载荷向量$ F_e^n $按“单元节点-全局节点”的自由度映射关系叠加至对应位置, 完成全局载荷向量$ F^n $的组装.其中$ F^n $的第$ k $个分量代表全局第$ k $个节点所对应的热源贡献.

对于二维无热源的非稳态热传导方程, 类似一维情形, 采用Crank-Nicolson格式进行时间离散, 得到如下全离散格式

$ \left( M + \frac{\alpha \Delta t}{2} K \right) u^{n+1} = \left( M - \frac{\alpha \Delta t}{2} K \right) u^n. $

由于二维问题在边界$ \partial\Omega $上施加了齐次Dirichlet边界条件即边界节点温度恒为0, 因此在求解时只需对内部节点进行未知量求解, 从而有效减少计算规模并提高效率.

为此, 我们从全局质量矩阵$ M $和全局刚度矩阵$ K $中提取对应于内部节点的子矩阵, 记为$ M_{\mathrm{int}} $和$ K_{\mathrm{int}} $.相应地, 定义内部节点的温度向量$ u_{\mathrm{int}}^n $和$ u_{\mathrm{int}}^{n+1} $, 分别表示$ t^n $和$ t^{n+1} $时刻所有内部节点上的温度值.由此得到仅关于内部节点的离散方程组

$ \begin{equation} A_{\mathrm{int}} u_{\mathrm{int}}^{n+1} = b_{\mathrm{int}}, \end{equation} $

(2.6)

其中

$ A_{\mathrm{int}} = M_{\mathrm{int}} + \frac{\alpha \Delta t}{2} K_{\mathrm{int}}, \quad b_{\mathrm{int}} = \left( M_{\mathrm{int}} - \frac{\alpha \Delta t}{2} K_{\mathrm{int}} \right) u_{\mathrm{int}}^n. $

各符号含义如下

● $ M_{\mathrm{int}}, \, K_{\mathrm{int}} $: 分别为全局质量矩阵$ M $和刚度矩阵$ K $中对应于内部节点的子矩阵;

● $ A_{\mathrm{int}} $: 内部节点系统的线性方程系数矩阵;

● $ b_{\mathrm{int}} $: 时间推进的右端项向量, 由$ t^n $时刻的内部节点温度及内部子矩阵计算得到;

● $ u_{\mathrm{int}}^n, \, u_{\mathrm{int}}^{n+1} $: 分别为$ t^n $和$ t^{n+1} $时刻的内部节点温度向量, 仅包含非边界节点的自由度;

● $ \alpha $: 热扩散系数;

● $ \Delta t $: 时间步长.

北太天元在本文中的作用并不局限于数值计算工具本身, 更体现在其对偏微分方程数值解实操教学全过程的支撑.其稀疏矩阵运算功能可高效处理有限元离散后的线性方程组, 科学可视化模块能够直观展示数值解、精确解及误差分布, 交互式编程环境便于课堂演示与算法调试, 相关PDE工具也有助于教师围绕典型模型快速组织实验案例.因此, 借助北太天元平台, 可以将离散推导、程序实现、结果展示与误差分析整合为较为连贯的教学过程, 从而为偏微分方程数值解课程的实操教学提供有效支撑.

首先展示北太天元自带的二维网格剖分功能及结果.如图 1所示, 展示的是二维求解域$ \Omega=[0, 1]\times[0, 1] $的剖分网格$ x $、$ y $方向节点数均取$ N_x=N_y=11 $, 三角形单元总数为$ N_e=2(N_x-1)(N_y-1) $, 网格步长满足$ h_x=h_y=\frac1{N_x-1}=\frac1{N_y-1} $.图 2为对应的网格剖分代码, 借助北太天元的内置工具箱可快速实现三角形网格剖分, 显著降低人工编程工作量.

一维问题结构相对简单, 便于学生首先理解有限元离散、时间推进与数值结果之间的对应关系, 因此可作为课堂实操教学的入门案例.根据(2.5)的离散格式进行一维求解, 求解区域$ \Omega=[0, 1] $, 单元数$ N_e=64 $, 网格长度$ h_e=\frac1{N_e} $, 时间步$ \Delta t=10^{-3} $, 求解$ t\in [0, 0.1] $的温度场.图 3是一维数值结果图, 其中子图(d)为$ t=0.1 $时的数值解与解析真解的对比, 可以看出数值解和精确解吻合.

在一维案例的基础上, 进一步考察二维问题, 有助于学生理解有限元方法在更复杂区域离散和更高维度计算中的实现特点, 从而增强其对算法普适性的认识. 根据(2.6)的离散格式进行二维求解, 求解区域$ \Omega=[0, 2]\times[0, 2] $.划分的网格中, $ x, y $方向节点数$ N_x=N_y=21 $, 总单元数$ N_e=2(N_x-1)(N_y-1) $, 网格步长$ h_x=h_y=\frac{2}{N_x-1}=\frac{2}{N_y-1} $, 时间步$ \Delta t=10^{-3} $, 求解$ t\in[0, 0.1] $的温度场.图 4为$ t=0.1 $时数值解、真解及绝对误差的三维可视化.子图(a)代表数值解, 子图(b)代表精确解, 从子图(c)可以看出, 最大绝对误差为$ 6.1865\times10^{-4} $.最后, 通过子图(d)可以看出区域中心点温度随着时间的推进, 精确解和数值解也吻合得很好.

在有限元方法的误差分析中, 收敛阶是衡量数值解逼近精确解速度的重要指标.对于线性有限元空间离散, 若真解$ u $充分光滑(例如$ u \in H^2(\Omega) \cap H_0^1(\Omega) $), 则其在$ L^2(\Omega) $范数下的误差满足如下最优阶估计：

$ \begin{equation} \| u - u_h \|_{L^2(\Omega)} \leq C h^2 |u|_{H^2(\Omega)}, \end{equation} $

(3.1)

其中$ u_h $为有限元近似解, $ h $表示网格的最大单元直径, $ C>0 $为与$ h $无关的常数, $ |u|_{H^2(\Omega)} $为解的二阶Sobolev半范数.该不等式表明, 在理想条件下, 线性有限元方法在$ L^2 $范数下具有二阶收敛性, 即收敛阶$ p = 2 $^[9].

本文所求解的二维热传导方程属于典型的抛物型偏微分方程, 其空间离散采用线性有限元法.理论上, 当时间步长足够小以抑制时间截断误差的影响时, 整体误差主要由空间离散主导, 因此应能观察到接近二阶的空间收敛行为.

为验证算法的空间精度, 选取计算区域$ \Omega = [0, 2] \times [0, 2] $, 并采用四组逐渐加密的结构化网格：$ 10\times10 $、$ 20\times20 $、$ 40\times40 $和$ 80\times80 $, 对应的最大网格尺寸分别为$ h = 0.2000 $、$ 0.1000 $、$ 0.0500 $、$ 0.0250 $.时间区间取为$ t \in [0, 0.1] $, 时间步长固定为$ \Delta t = 10^{-3} $, 以确保时间方向的离散误差远小于空间误差, 从而避免对空间收敛阶的干扰.

在各网格下求解后, 计算最终时刻的$ L^2 $范数误差, 并利用相邻网格对之间的误差比值估算实际收敛阶.二维收敛阶$ p_{\text{pair}} $的计算公式为

$ \begin{equation} p_{\text{pair}} = \frac{\ln\left(\frac{E_k}{E_{k+1}}\right)}{\ln\left(\frac{h_k}{h_{k+1}}\right)}, \end{equation} $

(3.2)

其中$ E_k $和$ h_k $分别表示第$ k $组网格下的误差与网格尺寸.

计算结果如表 1所示

结合表 1和图 5可知, 随着网格不断加密, $ L^2 $误差显著减小, 且逐对计算的收敛阶均非常接近2, 这一现象与线性有限元方法的理论预测高度一致.

从教学角度看, 收敛阶验证不仅用于说明算法正确, 也能够帮助学生建立“理论误差估计—数值实验现象—程序实现结果”三者之间的对应关系, 从而提高其对数值分析思想的整体理解.本文所采用的有限元格式在空间方向上达到了最优的二阶收敛精度, 验证了算法的正确性与有效性.该结果也说明, 在合理控制时间步长的前提下, 空间离散能够主导整体误差行为, 满足$ O(h^2) $的收敛速率.

本节给出一维与二维问题的关键代码模块.这里展示代码的目的, 并非单纯说明软件操作步骤, 而是用于呈现从网格生成、有限元矩阵组装到时间推进和误差分析的完整实现链条, 使学生能够将课堂中的离散公式与实际程序结构对应起来.换言之, 代码展示在本文中承担的是“连接理论推导与实验实现”的教学功能, 而不是简单的软件操作说明.其中, 一维问题的关键代码如图 6所示, 二维问题的关键代码如图 7所示.

本文以热传导方程有限元求解为案例开展教学实践, 其意义并不止于完成一个偏微分方程数值实验, 而在于借助北太天元平台重构偏微分方程数值解课程的实操教学组织方式.与传统以公式讲授和结果展示为主的教学模式相比, 本文更强调以具体问题为牵引, 将数学建模、离散推导、程序实现与误差分析整合为一个连续的学习过程, 使学生在“做中学”的过程中理解数值方法的形成机制与应用逻辑.基于这一思路, 北太天元在本文中发挥的作用并非单一的软件替代, 而是教学过程重构的技术支撑.

偏微分方程数值解课程具有较强的理论性与实践性.传统教学中, 教师通常需要在课堂上完成控制方程推导、离散格式说明和算法流程讲解, 而程序实现与结果验证往往被压缩到课后完成.对于初学者而言, 从弱形式到单元矩阵、从局部组装到全局求解、从数值解到误差分析, 这一链条较长且抽象, 学生容易停留在“知道结论”而难以真正完成“自主实现”.同时, 若实操教学依赖多个软件环境切换, 教师在课前准备、案例搭建和课堂演示上的投入也会明显增加, 进而影响课程实施的稳定性与持续性.

依托北太天元平台, 可以在较为统一的环境中完成网格生成、矩阵运算、程序编写与结果可视化, 从而将原本分散的教学环节整合起来.对于教师而言, 这种一体化环境减少了实操教学对额外软件配置和多平台切换的依赖, 降低了案例准备成本与课堂组织难度; 对于学生而言, 则缩短了从理论公式到程序结果之间的路径, 使其能够把更多精力放在理解有限元思想、把握离散步骤和分析数值结果上, 而不是消耗在复杂环境配置与机械性操作之中.由此, 课程实操环节更容易真正进入课堂教学, 而不再停留在课后附加任务的层面.

更为重要的是, 北太天元支撑下的实操过程能够形成清晰的任务分层.教师可先引导学生完成控制方程与边界条件的理解, 再组织其依次完成网格剖分、局部矩阵构造、时间推进求解及误差验证等任务, 使学生在逐步推进中建立对有限元方法的整体认识.这种分层式训练有助于缓解传统教学中“理论听懂了、程序不会写、结果看不懂”的典型问题, 提升课堂实操的落地性与教学效果.

将北太天元引入偏微分方程数值解教学, 其价值不仅体现在软件来源的国产化, 更体现在其能够支撑课程内容与实践训练的深度融合.传统教学中, 学生常将“数学理论”和“数值实验”视为彼此割裂的两个部分: 前者是课堂上的推导与证明, 后者则是课后的程序运行与图形展示.本文所构建的教学路径则尝试打破这种分离状态, 使数学模型、离散思想、代码表达与数值现象在同一教学过程中相互印证.

以热传导方程为例, 学生在学习过程中不仅需要理解控制方程、边界条件和初始条件的物理含义, 还需要进一步掌握有限元空间离散的基本思路、Crank–Nicolson格式的时间推进机制, 以及误差与收敛阶所反映的数值特征.当这些内容能够在同一平台中被连续呈现时, 抽象公式不再只是纸面上的推导结果, 而能够通过程序运行和图像输出来获得直观验证.这样一来, 学生对数值方法的认识将从“记住算法步骤”逐步转向“理解算法为何成立、为何有效、适用于何种问题”.

在这一过程中, 北太天元的意义在于为“理论—实现—验证”的闭环提供了操作基础.学生不仅能够看到离散格式如何转化为可求解的代数方程组, 还能够通过数值解与精确解的比较、误差曲线与收敛阶的分析, 进一步理解有限元方法的精度特征与适用范围.由此, 教学活动不再局限于数值实验任务的完成, 而是转向对学生数值建模能力、算法实现能力和结果分析能力的综合培养.这种融合式教学模式, 更符合当前高等教育对应用型与创新型人才培养的要求.

教学改革能否推广, 关键不只在于案例本身是否完整, 更在于实施路径是否足够清晰、门槛是否足够可控.本文围绕北太天元构建的教学方案具有较强的本土化特征, 其推广价值主要体现在两个方面.

其一, 该方案具有较好的课程适配性.热传导方程是偏微分方程数值解课程中的经典模型, 其所涉及的弱形式建立、有限元离散、时间推进与误差分析等内容, 具有较强的代表性.围绕这一案例形成的教学流程, 不仅适用于偏微分方程数值解课程, 也可迁移到数值计算方法、有限元方法、科学计算导论等相关课程之中.教师可根据课程层次和学时安排, 对案例规模、代码复杂度和实验要求作适当调整, 从而实现由基础实操到综合训练的分层应用.

其二, 该方案具有较强的实施可复制性.本文形成的教学路径并不依赖高度定制化的实验条件, 而是建立在较为明确的教学环节之上, 即问题提出、模型建立、离散推导、程序实现、结果验证与教学反思.对于多数高校而言, 只要具备基本的计算机实训条件, 便可以参照这一流程组织教学.这意味着北太天元不仅能够作为一个软件工具进入课堂, 更能够作为教学方案的一部分被整体引入, 实现从单次案例展示到稳定课程实践的转化.

从更广的意义上看, 国产科学计算软件进入教学场景, 不应只停留在“能替代”的层面, 更应进一步体现“好组织、好实施、好推广”的教育价值.本文所构建的教学路径, 正是试图在课程实践中回答这一问题: 即国产软件不仅可以承担数值计算任务, 也可以成为组织实操教学、提升教学效果、培养学生实践能力的重要支撑.

从教学目标看, 本文方案的价值主要体现在以下几个方面.首先, 它有助于提升学生对偏微分方程数值解全过程的理解深度.学生不再仅停留于对离散格式和程序代码的局部认识, 而是能够在完整案例中把握模型、算法与结果之间的内在联系.其次, 它有助于强化学生的实践能力.通过亲自完成模型离散、程序编写、数值求解与误差分析, 学生能够逐步形成将数学问题转化为计算问题的基本能力.再次, 它有助于解决课程教学中的实际难题, 特别是理论讲授与实操训练脱节、实验案例组织难、教学准备成本高等问题.最后, 它为同类课程提供了较为清晰的推广范式, 使国产软件能够以教学方案而非单纯工具的形式融入课堂.

因此, 本文的创新性并不单纯体现在完成了一维和二维热传导方程的数值实验, 而在于依托北太天元平台, 将偏微分方程数值解课程中的典型内容组织为一条可实施、可复制、可推广的教学路径.该路径既回应了国产软件融入课程教学的现实需求, 也为相关课程的实操教学改革提供了新的思路.

本文以一维和二维非稳态热传导方程为案例, 基于北太天元平台完成了有限元空间离散、Crank–Nicolson时间推进、程序实现与数值验证, 并在此基础上讨论了国产科学计算软件支撑偏微分方程数值解实操教学的实施路径.数值结果表明, 所构建的离散格式具有良好的精度与稳定性, 能够满足课程教学中的基本实验需求; 更重要的是, 相关教学实践表明, 借助北太天元可以较好地打通理论讲授、程序实现与结果分析之间的联系, 使学生在完整的问题求解过程中理解数值方法的核心思想.

与将软件仅作为实验工具的做法不同, 本文更强调北太天元在教学组织中的平台作用.依托其矩阵运算、可视化展示和程序开发等功能, 可以将偏微分方程数值解课程中的抽象内容转化为层次清晰、步骤明确的实操任务, 从而降低教师实施课程实验的准备成本, 增强学生参与计算实践的可达性与获得感.由此, 北太天元的应用价值不仅体现在完成数值实验, 更体现在支持课程教学从“结果展示”走向“过程训练”, 从“软件使用”走向“能力培养”.

同时, 本文形成的教学路径具有较强的本土化与可推广特征.围绕典型偏微分方程案例建立起来的“理论讲授—案例驱动—代码实践—结果验证”教学流程, 既适用于本课程的课堂教学, 也可为数值计算方法、有限元方法及相关实操课程提供借鉴.因而, 本文工作的意义不仅在于给出一个基于国产软件的教学案例, 更在于为国产科学计算平台融入高等数学类课程教学提供了一种可复制的实施方案.

后续研究可进一步结合不同层次学生的学习特点, 引入更多类型的偏微分方程模型与更丰富的教学评价方式, 持续完善基于国产软件平台的偏微分方程数值解教学体系, 进一步发挥国产工具在课程建设、实践教学与创新人才培养中的作用.