8.3 求解显式差分格式的方法
例8.1 一维非瞬态对流扩散方程（8.13）中，求解域为[0，1]，初始条件为T（x，0）=sin（2πx），左边界条件为T（0，t）=-
），右边界条件为T（1，t）=
（2π（1-αt）），方程（8.13）的解析解为
1.静态模板的制作
第一步，输入已知数据。在单元格A2、B2和C2中分别输入区间的长度、流速和扩散系数（这里分别输入“1”“2”和“0.1”）；在单元格D2、E2、F2、G2、H2和I2中分别输入空间步长、时间步长、常数λ=
、计算停止时间的公式和计算总计算次数的公式（这里分别输入“0.05”“0.001”“=B2*E2/D2”“=C2*E2/（D2^2）”“0.5”和“=H2/E2”），如图8.1所示。
图8.1 输入已知数据
第二步，计算解析解。在单元格B5中输入区间的左边界“=0”，在单元格C5中输入“=B5+$D$2”；在单元格B6中用式（8.15）输入解析解计算公式“=EXP（-4*$C$2*PI（）^2*$H$2）*SIN（2*PI（）*（B5-$B$2*$H$2））”；选中单元格B6，把鼠标指针移动到单元格B6的右下角，当其变成“
”符号时，按住鼠标左键向右拖动一列；选中单元格区域B5：B6，把鼠标指针移动到单元格B6的右下角，当其变成“
”符号时，按住鼠标左键向右拖动到单元格V6，至此就完成了在区间[0，1]内，每一个节点处解析解的计算，如图8.2所示。
图8.2 解析解
第三步，计算数值解。在A9中输入“0”；在单元格区域B9：V9中输入初始条件，在单元格B9中输入“=SIN（2*PI（）*（B5））”，然后向右拖动到单元格V9即可。在单元格A10中输入“=A9+1”，在单元格B10中输入左边界计算公式“=EXP（-4*$C$2*PI（）^2*（$E$2*A10））*SIN（2*PI（）*（$B$5-$B$2*（$E$2*A10）））”，在单元格C10中用式（8.14）输入“=C9+$G$2*（B9-2*C9+D9）-0.5*$F$2*（D9-B9）”，然后向右拖动到单元格V10，就得到第一步的数值解。选中单元格区域A10：V10，把鼠标指针移动到单元格V10的右下角，当其变成“
”符号时按住鼠标左键向下拖动一列，就得到下一步的数值解，要计算50步的话需按住鼠标左键拖动到单元格V59，得到t=0.05时刻的数值解。
第四步，绘制图像。选中单元格区域B5：V6，依次单击“插入”、图表下拉按钮、“更多散点图
）...”，选择“散点图→带平滑线的散点图”，单击“确定”按钮，得到解析解的图像；选中解析解的图像，单击右键，选择“选择数据”，打开“选择数据源”对话框，单击“添加”按钮，弹出“编辑数据系列”对话框后，在系列名称内填“数值解”，在X轴系列值内选择填“=Sheet2！$B$5：$V$5”，在Y轴系列值内选择填“=Sheet2！$B$59：$V$59”，单击“确定”按钮且对图进行装饰，得到图8.3所示的比较图。
图8.3 解析解与数值解的比较图
2.动态模板的制作(https://www.daowen.com)
第一步，设计迭代计算。在“文件”菜单下选择“选项”，在“选项”对话框中选择“公式”副对话框。选择启用迭代计算，最多迭代计算次数设为1，最大误差设为0.001，单击“确定”按钮。
第二步，输入已知数据。在单元格A2、B2和C2中分别输入区间的长度、流速和扩散系数（在这里输入“1”“2”和“0.1”）；在单元格D2、E2、F2、G2、H2和I2中分别输入空间步长、时间步长、常数λ=
、计算停止时间的公式和计算总计算次数的公式（在这里输入“0.05”“0.001”“=B2*E2/D2”“=C2*E2/（D2^2）”“1”和“=H2/E2”）；在空单元格J2中输入实现循环的表达式“=IF（J2＜I2/2，J2+1，1）”；在单元格K2中输入计算正在计算时刻的表达式“=2*J2*E2”，如图8.4所示。
图8.4 参数和循环迭代控制语句的输入
第三步，输入初始条件和边界条件。在单元格B5中输入区间的左边界“=0”，在单元格C5中输入“=B5+$D$2”；在单元格B6中用初始条件表达式输入“=SIN（2*PI（）*（B5））”。选中单元格B6，把鼠标指针移动到单元格B6的右下角，当其变成“
”符号时，按住鼠标左键向右拖动一列，再选中单元格区域B5：B6，把鼠标指针移动到单元格B6的右下角，当其变成“
”符号时，按住鼠标左键向右拖动到单元格V6，完成在区间[0，1]内每一个节点初始值的计算。在单元格A6、A7和A8中分别输入显示计算步数公式“=＂k=0＂”“=＂k-1=＂&2*J2-1”和“=＂k=＂&2*J2”；用左边界条件表达式在单元格B7和B8中分别输入“=EXP（-4*$C$2*PI（）^2*（K2-E2））*SIN（2*PI（）*（B5-$B$2*（K2-E2）））”和“=EXP（-4*$C$2*PI（）^2*K2）*SIN（2*PI（）*（B5-$B$2*K2））”；用右边界条件表达式在单元格V7和V8中分别输入“=EXP（-4*$C$2*PI（）^2*（K2-E2））*SIN（2*PI（）*（V5-$B$2*（K2-E2）））”和“=EXP（-4*$C$2*PI（）^2*K2）*SIN（2*PI（）*（V5-$B$2*K2））”。
第四步，计算第一步和第二步的数值解。在单元格C7中用式（8.14）输入“=C6+$G$2*（B6-2*C6+D6）-0.5*$F$2*（D6-B6）”。选中单元格C7，把鼠标指针移动到单元格C7的右下角，当其变成“
”符号时，按住鼠标左键向右拖动到单元格U7，就得到第一步的数值解。选中单元格区域C7：U7，把鼠标指针移动到单元格U7的右下角，当其变成“
”符号时，按住鼠标左键向下拖动一列，就得到下一步的数值解。
第五步，实现循环计算。为了实现循环计算，将单元格C7中的计算公式改为“=IF（$J$2=1，C6+$G$2*（B6-2*C6+D6）-0.5*$F$2*（D6-B6），C8+$G$2*（B8-2*C8+D8）-0.5*$F$2*（D8-B8））”。选中单元格C7，把鼠标指针移动到单元格C7的右下角，当其变成“
”符号时，按住鼠标左键向右拖动到单元格U7。这时，如果我们按功能键F9，则在单元格J2的值每一次增加1，循环在1～500之间，当J2的值不等于1时，单元格区域C7：U7用单元格区域C8：U8的值进行计算，而单元格区域C8：U8的值随着单元格区域C7：U7值的更新而更新，这样在单元格区域B8：V8内就得到单元格K2的值在对应时刻的数值解，如图8.5所示。
图8.5 动态数值解
第六步，计算解析解。在单元格区域B11：V11中用式（8.15）输入计算解析解公式，选中单元格区域B12：V12，输入解析解与数值解绝对误差的计算公式“=ABS（B11：V11-B8：V8）”，同时按住“Ctrl+Shift+Enter”键并确认。
第七步，绘制图像。在一个坐标系内同时绘制解析解与数值解图像，在另一个坐标系内绘制误差图。图8.6（a）是当t=0.2（即第k=100）时解析解与数值解的比较图，图8.6（b）是误差图。
图8.6 当t=0.2时解析解与数值解的比较图和误差图
