抛物型方程差分格式实验

状态：由《Matlab求解数值结果.pdf》12 页及 Parabolic_1D.m、Paradomic_2D.m 整理。PDF 中的结果表已复核；本轮没有 MATLAB/Octave 环境，源码只完成静态检查，未重新执行。

1. 一维热方程：稳定性比格式名称更重要

考虑

\[ u_t=u_{xx}+\sin t, \quad 0< x<1, \qquad u_x(t,0)=u_x(t,1)=0, \]

初值与精确解为

\[ u(x,0)=\cos(\pi x), \qquad u(x,t)=e^{-\pi^2t}\cos(\pi x)+1-\cos t. \]

令

\[ r=\frac{\tau}{h^2}. \]

实验比较三类格式：

PDF 记录的 \(t=1\) 采样结果显示：

格式与网格	\(r\)	采样点最大绝对误差	现象
显式，\(h=1/40,\tau=1/3200\)	\(0.5\)	约 \(1.32\times10^{-4}\)	稳定
显式，\(h=1/80,\tau=1/12800\)	\(0.5\)	约 \(3.3\times10^{-5}\)	网格减半后误差约缩小四倍
显式，\(h=1/80,\tau=1/3200\)	\(2\)	第 40 步已到 \(10^{17}\) 量级	违反稳定条件后爆炸
隐式，\(h=1/40,\tau=1/1600\)	\(1\)	约 \(2.65\times10^{-4}\)	稳定但更耗散
隐式，\(h=1/80,\tau=1/3200\)	\(2\)	约 \(1.32\times10^{-4}\)	大 \(r\) 下仍稳定
Crank-Nicolson，两组网格	\(1,2\)	六位小数下显示为 \(0\)	记录精度内与真解非常接近

“显示为零”只说明误差小于表格的六位小数分辨率，不能据此声称误差严格为零。

第二个实验考虑

\[ u_t=\frac1{16}(u_{xx}+u_{yy}), \qquad (x,y)\in(0,1)^2, \]

配合 \(x\) 方向齐次 Dirichlet 边界、\(y\) 方向复制型离散 Neumann 边界，以及初值

\[ u(x,y,0)=\sin(\pi x)\cos(\pi y). \]

精确解为

\[ u(x,y,t)=\sin(\pi x)\cos(\pi y)e^{-\pi^2t/8}. \]

实验比较：

在 \(h=1/40,\tau=1/1600,t=1\) 的九个采样点上，PDF 记录的最大绝对误差为：

这些单网格结果只能比较当前设置下的数值表现。源码还会在 \(N=20,40,80\) 上计算无穷范数误差和经验收敛阶，但 PDF 没有保存最终阶数表，因此这里不补写未实际看到的数值。

二维隐式格式若直接求解，会形成 Kronecker 和结构的大型稀疏系统。ADI 与 LOD 利用坐标方向可分性，把二维问题化成多组一维带状系统：

这正是工业计算页面中“离散结构决定求解器与硬件映射”的具体例子。

脚本会输出 CSV、数值/精确解对照和误差图。二维源码的本地输出目录名已从原文件的 Paradomic_2D_output 拼写修正为 Parabolic_2D_output；数值算法没有改写。