8.7 求解隐式差分格式的方法
例8.3 一维扩散方程公式（8.23）中，求解域为[0，1]，初始条件为
左边界条件和右边界条件为
它的解析解为
就本问题来说方程（8.21）中的bn为零，不必计算，此时方程（8.29）变为
下面介绍动态模板的制作过程。
第一步，启动迭代计算。在“文件”菜单下选择“选项”，在“选项”对话框中单击“公式”副对话框。选择启用迭代计算，最多迭代计算次数设为1，最大误差设为0.001并单击“确定”按钮。
第二步，输入常数项和相关参数。按照图8.11，在打开的Excel表格中分别输入常数α，求解域区间[a，b]，加权系数θ，x方向的网格步长h，λ和时间步长τ的值。当输入时间步长时，在单元格G2中按差分格式的稳定性条件输入“=IF（D2＜0.5，E2^2/（2*A2）*1/（1-2*D2），E2^2*F2）”。在单元格H2中输入计算终止时间T；在单元格I2中输入计算总次数n的公式“=H2/G2”，在单元格J2中输入循环迭代控制值的语句“=IF（J2＜I2/2，J2+1，1）”；在单元格K2中输入计算现在时刻的语句“=2*J2*G2”。
图8.11 参数和循环迭代控制语句的输入
第三步，输入矩阵A和矩阵B。用矩阵A的计算公式在单元格A5中输入“=1+2*A2*F2*D2”；在单元格B5中输入“=-A2*F2*D2”；在单元格A6中输入“=B5”；选中单元格区域B6：C6，输入“=A5：B5”，同时按“Shift+Ctrl+Enter”键确认；选中单元格区域B7：D7，输入“=A6：C6”，同时按“Shift+Ctrl+Enter”键确认；用类似方法输入矩阵A的其他元素，得到图8.12所示的结果。同理，用矩阵B的计算公式输入矩阵B的元素，得到图8.13所示的结果。
(https://www.daowen.com)
图8.12 矩阵A
图8.13 矩阵B
第四步，输入初始条件。在单元格A27中输入“=B2”，在单元格A28中输入“=A27+$E$2”。选中单元格A28，把鼠标指针移动到单元格A28的右下角，待其变成“
”符号时，按住鼠标左键向下拖动到单元格A的值变到左边界b的值1为止，松开鼠标左键。在单元格B27中输入“=SIN（PI（）*A27）”，用上述同样的选择和拖动方法计算每一点的初始条件值，如图8.14所示。
图8.14 数值解与解析解的比较
第五步，输入边界条件。单元格C27、D27、C37和D27中相应的边界条件或边界条件的计算公式在本例中全为零，如图8.14所示。
第六步，输入第一步和第二步的公式。选中单元格区域C28：C36，用式（8.23）输入计算第一步的公式“=MMULT（MINVERSE（K5：S13），MMULT（A16：K24，B27：B37））”，同时按“Shift+Ctrl+Enter”键确认；选中单元格区域D28：D36，输入计算第二步的公式“=MMULT（MINVERSE（K5：S13），MMULT（A16：K24，D27：D37））”，同时按“Shift+Ctrl+Enter”键确认，如图8.14所示。
第七步，输入循环迭代计算语句。选中单元格区域C28：C36，输入迭代计算公式“=IF（M2=1，MMULT（MINVERSE（K5：S13），MMULT（A16：K24，B27：B37）），MMULT（MINVERSE（K5：S13），MMULT（A16：K24，E27：E37）））”，同时按“Shift+Ctrl+Enter”键确认，如图8.14所示。
第八步，计算解析解和误差。用式（8.30）在单元格E27中输入“=SIN（PI（）*A27）*EXP（-POWER（$A$2*PI（），2）*$K$2）”，用上述同样的方法计算相应时刻的解析解；选中单元格区域F27：F37，输入绝对误差计算式“=ABS（D27：D37-E27：E37）”，同时按“Shift+Ctrl+Enter”键确认，如图8.14所示。
第九步，绘制图像。在一个坐标系内同时绘制解析解与数值解图像，在另一个坐标系内绘制误差图，图8.15（a）是t=0.2（即第k=100）时解析解与数值解的比较图，图8.15（b）是误差图。
图8.15 动态图像的比较图和误差图
以上操作完成以后，按住功能键F9，随着单元格J2值的循环变化，相应时刻的数值解和解析解也变化。如果改变加权系数θ，可以观察加权系数θ对数值解的影响。
