永冻土层热传导问题数学建模

1、摘要本文针对永冻土层上关于路基热传导的问题，通过对不同材料层的密度、比 热容、传热系数进行研究，建立微分方程模型，利用 Matlab与Lingo软件进行 求解。问题一，考虑在分析各材料层进行后，给出空气温度传入路基规律，以及各 材料层的温度分布。同时，已知外界的温度是关于时间的函数， 冻土层的温度是 不变的零下温度。首先，我们通过中国选矿技术网以及中国天气网分别获得各材 料层的密度、比热容、传热系数等数据，和拉萨最近 24小时的温度数据。通过 拟合得到温度与时间的关系函数，建立一维热传导方程的微分方程模型。 随后利 用向前差分的方法求出方程的近似数值解， 因为界面处的热传导率处于平衡，且 温度

2、相等，则可以一层一层向下计算得出各材料层的温度分布规律。问题二，考虑在一些设备的支架不能固定在解冻土层上， 必须固定在永冻土 层中的情况下，地下土层的解冻位置，并给出解冻砂土与冻结砂土的分界线。 由 问题一的求解可以计算出L4的值，在已知上界x值与下界的0C温度值后，同问 题一的求解方法可以给出解冻砂土与冻结砂土的分界线。问题三，考虑结合温度分布、成本及耐用性，给出各层材料的最佳厚度。结 合铁路建设施工保障，我们将耐用性作为出发点，分别从压强及压实度考虑耐用 性的约束条件。由于压实度与含水量存在联系，而含水量与温度存在关系，故建 立起压实度与温度的相关关系。将成本作为目标函数，压实度与压强的限

3、制作为 约束条件，建立线性规划模型，由此解出最少成本为 46334.72元。问题四，考虑在以上问题的基础上，结合我国青藏铁路永冻土层地基进行仿 真，并为施工单位提出合理建议。因为本题的前三问即是在查阅青藏铁路路基修 建相关数据的基础上进行的，故问题四的仿真即已经得到相应的解决。 通过对以 上问题的求解进行合理性分析，即可对施工单位给出合理建议。为了简化计算量，提高求解速度，本题中的微分方程模型使用向前差分的方 法求出近似数值解，而且对模型的可行性及有效性进行了一定的分析， 所得结果 十分合理。本文的优点在于利用差分方法求解一维热传导微分方程模型的近似数值解， 使得材料内界面的条件处理得较为容易

4、。 同时，在前三个问题的求解中，将问题 背景设定在青藏铁路的修建中，在一定程度上对问题四的求解提供了较大的帮助。关键词：热传导问题 抛物型方程 数值模拟1问题重述在永冻土地上铺设道路、飞机跑道和某些结构的地基。分析这类地基的结构, 有沥青层、混凝土层、干砂层、石子层、绝缘材料层，再下面是湿的沙土层和冻 土层（参见图1）。已知，外界的温度是关于时间的函数，冻土层的温度是不变 的零下温度。请回答以下问题：沥青层混凝土层干砂石层防水隔温层半冻砂土层图1:永冻土地地基结构图（1）分析空气温度传入路基规律，各材料层的温度分布；（2）由于一些设备的支架不能固定在解冻土层上， 必须固定在永冻土层中, 因此确

5、定地下土层的解冻位置非常重要，请给出解冻砂土与冻结砂土的分界线；（3） 请给出各层材料的最佳厚度（结合温度分布、成本和耐用性）；（4）请结合我国青藏铁路永冻土层地基进行仿真，并给施工单位给出合理 建议。2问题分析本问题要求建立模型和设计算法，在永冻土地上铺设道路、飞机跑道和某些 结构的地基的问题背景中，结合温度分布、成本和耐用性，给出各层材料的最佳 厚度，并对我国青藏铁路永冻土层地基进行仿真。 由于问题的局限性，我们首先 应该找到合理的数据来对相关背景进行一定的了解， 定位模型的求解方法，循序 渐进地对此问题进行合理的求解。问题一中，要求在题中所给的情况下，分析空气温度传入路基规律，以及各 材

6、料层的温度分布。经过分析，我们发现建立模型所需的数据不足， 故应首先搜 集题目中各材料热传导相关的密度、导热系数、比热容以及环境温度等数据。问题二中，要求确定地下土层的解冻位置，给出解冻砂土与冻结砂土的分界 线，以应对由于一些设备的支架不能固定在解冻土层上，必须固定在永冻土层中 的情况。由问题一建立的模型，本问题即是求出温度为 0 C时所对应的距离。问题三中，要求结合温度分布、成本和耐用性，给出各层材料的最佳厚度。 在分析问题一、二后，则是对各目标进行权重约束，以求得最佳的厚度。问题四中，要求结合我国青藏铁路永冻土层地基进行仿真， 并给施工单位给 出合理建议。3模型假设1、假设环境温度不会发生

7、突变，没有极端天气出现；2、假设路基材料分布均匀；3、假设考虑压实度时，可以将半冻砂土层以上的四层材料合并为一层;4、假设对于各层只考虑重力所产生的压强。4符号说明a(t)外界的温度关于时间的函数To冻土层的温度K(i =1,2，,5)每个材料层的热传导率q每个材料层的比热R每个材料层的密度u要求解的第i个材料层的温度分布x地下某一点距地面的距离t时刻值li每层材料的厚度其他局部符号在引用时给出了具体说明。5模型的建立与求解5.1模型一的建立与求解问题一中要求分析空气温度传入路基规律，以及各材料层的温度分布。我们知道，一切稳定的数学物理问题都可以用椭圆型微分方程模型来描述，若在稳态之前有一个过

8、程，则讨论这种渐进稳定的数学物理过程可以用抛物型为方程模型 来描述1 05.1.1模型一的建立分析题目可知，地基的结构有沥青层，混凝土层，干砂石层，防水隔温层，半冻砂土层，外界的温度是关于时间的函数 a(t)，冻土层的温度是不变的零下温 度T0，每个材料层都有热传导率ki(i =1,2,5)，问题一则要分析空气温度传入路基的规律，以及各材料层的温度分布2。首先，由于路基各材料层是均匀的，所以要分析的热传导方程可归为一维的 热传导方程研究。路基各材料层则化为一维热传导方程为:u:ti2(1)K qCi'i；：x2其中，Ui是要求解的第i个材料层的温度分布；ki，q，一为此材料层的传热系数

9、，比热和密度，它们都是已知的常数，且上.0。Ci R由模型假设可知，层与层的界面处传导率处于平衡，且温度相等，即ki -= ki 1 U,Ui = Ui i, x = Li,（i = 1,2,3,4） （2）x:X基于以上分析，在半冻砂土层假设为湿砂层（k5是其传热系数）的情况下，所建立的路基热传导微分方程模型为* 2cui ki c uifit Ci R ex2ki 土 = ki Ul-l ,Ui 二 Ui i,x = Li,（i =1,2,3, 4）（3）exexUo（O,t） =a（t）iU5（ L5,t） =T05.1.2模型一的求解由以上所建立的抛物型微分方程模型进行分析，对问题一的

10、求解我们给出了如下的具体步骤：（1）依据题意，我们在中国选矿技术网中找到沥青层，混凝土层，干砂石层，防水隔温层和半冻砂土层的密度，比热容和热导系数等值，部分难以确定的材料选取了最接近的材料进行替代，所找到的数据材料如下表。表1各层材料相关参数表各材料层干密度吒（Kg /m3）热导系数九2（w/ m 旳比热容C(KJ/(Kg*)价格P(元 /t)沥青层（石油沥青）14000.271.682800混凝土层（c30混凝土）23001.510.92770.5干砂石16000.581.01705防水隔温层（乳化沥4000.121.551000青膨胀珍珠岩）砂土层（湿砂土）23501.421.07500以

11、上则为各材料层的密度、热导系数、比热容和价格（价格因素在问题三中 需要考虑，故在这里一起给出）的数据。此外，由于温度的传导涉及到冻土地带的环境温度，故我们结合青藏铁路的修建，选取拉萨最近 24小时内的温度3 考虑到中国天气网只可以给出最近 24小时的较精确温度，故我们问题的季节背 景选择在夏秋之交4。表2拉萨24小时各时刻温度表时刻000102030405060708091011温度121212131313131313121313时刻121314151617181920212223温度151718192122222120161312以上则为拉萨最近24小时各时刻的温度表。（2） 通过查阅铁路建

12、设国家规范，我们将沥青层的厚度设定为4cm,混凝土层 设定为10cm,干砂石层设定为7cm,防水隔温层设定为10cm5。（3）通过拟合数据得到外界温度关于时间的关系式，根据问题一中建立的第一 种边界的一维抛物型方程，考虑到解析解不易求得，故在这里我们利用差分的方 法求出其近似的数值解。差分求解抛物型偏微分方程的过程则是首先对 xt平面进行网格剖分，分别取h，.为x方向和t方向的步长，用两族平行直线x二Xk二kh（k = _1,_2,），t二tj二j （j =0,1,2,），将xt平面剖分成矩形网格，节点为（Xk,tJ，为简单起见， 记（k, j） =（Xk，tj），u（k, j） =u（x&l

13、t;,tj），则在网格内点（k,j）处，对一采取向前差商ct-2公式，对弓采用二阶中心差商公式，可得到一维热传导方程的差分近似：u(k, j 1)-u(k, j)k u(k 1) -2u(k) u(k -1)叮ih2-0（4）（4）求出第一层底部的温度后，令k =1，依据公式（2）即可求得下一层底部的 温度，以下每一层材料温度的求解依次循环以上步骤。5.1.3模型一的结果及分析针对5.1.2中的求解步骤，运用 Matlab编程求解。针对附录一（1）中的程 序文件，运行可得到如下结果。外界环境温度随时间的变化函数为：a(t) = 15.11 + 2.216cos(1.757t) + 3.687s

14、in(1.757t) -1.895cos(3.514t) + 1.728sin(3.514t)外界温度随时间变化的函数如下图:(5)6图1外界温度随时间变化的函数第一层求解所得的温度随距离或时间变化的图像分别为:图2第一层温度随距离分布图图3第一层温度随时间分布图其他层温度随时间，随距离的变化图像见附录一（2）。由上可知，当固定另8一因素时，温度随着与路基表面距离的增加越来越小， 其随时间的变化则呈现先 增加后减少的趋势。5.2模型二的建立与求解根据题意，问题二中要求确定地下土层的解冻位置，给出解冻砂土与冻结砂 土的分界线，以应对由于一些设备的支架不能固定在解冻土层上，必须固定在永 冻土层中的

15、情况。由问题一建立的模型，本问题即是求出 0 C时所对应的距离。5.2.1问题二模型的建立根据问题一中所建立的模型，即可得出第五层温度随着距离变化的图像关系。 由图分析，即可给出0 C时所对应的距离7。5.2.2问题二模型的求解问题二中建立的模型用来确定第五层的解冻位置，即在已知上界的 x值与下 界温度后，给出解冻砂土与冻结砂土的分界线。 由函数图象分析，即可给出0 C 时所对应的距离。5.2.3问题二模型的结果分析针对5.2.2中的求解步骤，依据问题一中所得的第四层温度随距离的图像分 析进行求解。所得图像如下图：:孟宾芬韦珂強E寸同变鱼t空间变逐H图4第五层温度随距离变化图像解得解冻砂土与冻

16、结砂土之间的分界线为路面下120cm-140cmo结合国家公路施工技术规范，分析数据可知所得结果比较合理，故可作为问题而的求解答案。5.2.4模型二的误差分析因为差分法是对数值解得近似，所以对最后结果的误差分析十分重要。 依据 附录二中的程序文件，得出如下的误差分析图像：埶？in与柚調弹比厂堆寫廉-0瞞03?丄-025乓070 15D.1DOS0J ai 0.20.3-0.40.5D.&07 D.B D.91x图5数值解与精确解比较图像由上可知，数值解与精确解的最大相对误差非常小，为0.42877%，且数值解的求解结果是稳定收敛的，故该结果可以作为本问题的求解。5.3模型三的建立与求解

17、由题目可知，问题三需要在问题一、二的求解启发中，结合温度分布、成本 和耐用性，给出各层材料的最佳厚度。在分析问题一、二后，则是对各目标进行 权重约束，以求得最佳的厚度。5.3.1问题三模型的建立由于题目中所给相关数据较少，由此我们思考对最佳厚度设置约束条件。 对 于耐用性，我们从压实度和压强两方面进行考虑。首先将上面四层材料合为一层 以简化计算，查阅文献可知：压实度=试验Hr100%'试验干密度二实际密度（1+含水率）而实际干密度则与温度、水分等条件相关，则在此时将温度分布考虑进模型 的建立。由此，需要拟合温度与水分的关系8。而对于压强，我们假设这里路基的承重来源只有上部材料的重力，则

18、可得出如下推导公式：p = F =沁=已也（ 6）S SSS对上述公式中出现的局部符号说明如下:P表示某层材料所承受的压强； F表示某层材料所承受的压力；S表示承压接触面的表面积；V 表示表面积 S对应的材料的体积; m表示体积 V对应的材料的质量；'表示材料的密度；l表示材料层的厚度。将目标函数设定为压实度大于 95%且能承受最大压强为2760MPa时求得的最小成本。用Q表示成本,min Q=1400 2.8h+2300 16712+1600 0.0713+400 II4 （7）同时添加对各层压强的限制，使各层所承受压强大于路基建设所应承载的最 小压强值。同时，分析数据知，1压实度=

19、（含水率）100% 95%（ 8）基于以上分析，用r表示含水率，以（7）为目标函数，（8）为约束条件建立 优化模型如下：min Q=14002.81+2300>M67l2+1600x0.07l3+400>dl4(9)1400x101, ±700 2300疋1012 启 2760 s.t <16001013 生 960400域1014 3 4001>0.95+r5.3.2问题三模型的求解由于问题三是在结合温度分布、成本、耐用性多种因素的基础上对厚度问题 进行的研究，所以需要考虑各个因素的相关数据，以及各个因素之间的相互关系， 具体的求解步骤如下：（1）查阅资料对

20、温度与水分含量的关系给出合理的拟合关系，并得出具体 函数公式；（2）对压强和压实度的影响因素进行化简，得到相应的约束关系；（3）列出压实度大于95%时的目标函数，在Lingo中求解线性规划问题即 可。5.3.3问题三模型的结果分析针对5.3.1中的模型，按照5.3.2中的求解步骤，首先利用Matlab求解得出温度随水分变化的函数关系式为：y =-4.139+33.95cos(0.01754x)+3.358si n(0.01754x)( 10)所画出的函数图像为:图6温度与水分的拟合图像拟合度数据如下：gof =sse: 0.6035rsquare: 0.9991dfe: 5adjrsquare

21、: 0.9985rmse: 0.3474有数据可以看出，温度与水分的函数关系拟合得较为合理。对所列出的目标函数与约束条件，运用Lingo编程求解，解得压实度大于95% 时的目标函数所取得最小值为46334.72元。此时，第一、二、三、四层的厚度分别为5cm，12cm，6cm, 10cm。同上分析，可以看出此数据与实际施工情况中 的数据相差较小，故可作为本题的求解。5.4模型四的建立与求解由题目可知，问题四要求结合我国青藏铁路永冻土层地基进行仿真，并给施工单位给出合理建议。在对之前问题求解的基础上，得出合理的分析与建议。5.4.1问题四模型的建立问题一、二中选取的数据即是青藏铁路路基建设的相关数

22、据，同时将外界环 境温度设定为拉萨市最近24小时的温度。故由此所得的厚度范围即为对青藏铁 路永冻土层地基建设的仿真9。5.4.2问题四模型的求解由上所述原理对问题四进行求解，具体求解步骤如下：（1）由以上三个问题的求解，给出各材料层最佳厚度取值表;（2）通过对所得结果进行分析，给出施工工人的合理化建议。 5.4.3问题四模型的结果分析由以上问题的求解，可以得出如下结果：表3各材料层最佳厚度取值表材料层厚度沥青层（石油沥青）5cm混凝土层（C30混凝土）12cm干砂石6cm防水隔温层（乳化沥青膨胀珍珠岩）10cm通过对数据与国家青藏铁路路基建设进行比对， 以上所得求解数据基本符合 铁路建设规范，

23、故可以作为本问题的仿真结果。543问题四模型求解对施工工人的建议（1）由以上求解可以给出各层材料大致的厚度取值，可将此作为一个基本 的施工参照给施工工人一定的理论指导；（2）由于模型建设中提出了相应的假设，故在实际问题的施工时，相关施工单位可结合自身施工经验对厚度取值进行一定的调整，同时也应进行多次项目预试验，以充分考虑铁路承重、修建成本以及后发事故。6模型的评价与改进方向6.1模型的评价由于本问题中，模型建立的背景为青藏铁路的修建，同时对于各材料的相关 热学参数进行了唯一取值，对于此情况下的铁路路基热传导问题可以合理有效地 解决。如果修建环境发生改变便，则各材料层成分以及环境温度均会发生改变

24、， 故应对该模型进行改进继而求得新的求解情况。6.1.1模型的优点（1）本模型中，所选取的数据均为与各材料层成分接近程度最大的材料的 热学参数，建立的线性模型函数设置合理，对于题中所给问题可以得到满意的求 解方案。（2）算法设计较为合适，计算量较小，能够高效运行出结果。（3）利用差分代替微分的思维求解一维热传导微分方程模型，可以使材料 内界面的条件处理得较为容易。（4）本模型在前三个问题的求解中，将问题背景设定在青藏铁路的修建中， 在一定程度上对问题四的求解提供了较大的帮助。6.1.2模型的缺点（1）问题一中，我们所选取的材料数据并不能充分考虑各层材料成分，温 度数据也仅为一个地点24小时内的

25、数据；（2）问题三中，在考虑压实度时，将半冻砂土层以上材料压缩为一种材料， 在一定程度上简化了问题模型。6.2模型的改进方向（1）对于所选取的温度及材料热学数据，可以选取多组求平均值，以使求 解结果更具代表性；（2）求解一维热传导微分方程时，对于初始条件的选取稍有欠缺，可以结 合实际情况选择更好的初值条件。对于差分求解一维热传导微分方程，虽在一定程度上减少了求解困难， 但其步长选取不太合理，可以进一步研究问题模型，以选取更合理地步长，以使 求解结果与真实值更加接近。14参考文献1 张辰熙，于明季节冻土层下人工冻土墙温度场实验研究J.低温建筑术.2010(6):144-146.2 文斌,吴青柏,

26、蒋观利,张鹏.模拟退火优化算法的冻土热传导参数反分析J.岩 土力学.2013(8):15-16.3 王建伟，李书民，李振卿.冻土层等效导热系数的探讨J.防渗技 术.2001 (3):36-38.4 李述训,吴通化.冻土温度状况研究方法和应用分析J.冰川冻土 .2004(4):142-148. 候彬彬，卢德唐，陈永辉,孔祥言.非达西渗流模型模拟青藏铁路抛石路基冻 土层温度J.水工渗流研究与应用进展第五届全国水利工程渗透学术研讨 会论文集.2006.史册.热传导方程有限差分法的MATLAB实现J.咸阳师范学院学 报.2009:24-26.7 张世民.青藏铁路多年冻土路基热-力稳定性数值仿真分析J.

27、 土木建筑与 环境工程.2012(34):3-7.8 穆彦虎，马巍,牛富俊,李国玉,王大雁，刘永智.青藏铁路多年冻土区普通路 基热状况监测分析J.冰川冻土 .2014(2):10-13.田立慧，凌贤长，王力娜,张锋.青藏铁路高温多年冻土区列车行驶路基长期 永久变形数值模拟研究J.地震工程学报.2014(4):21-25.15附录问题一 外界环境温度随时间变化的函数关系式拟合程序creatFit.m function fitresult, gof = createFit(x, y) %CREATEFIT(X,Y) % Create a fit.% Data for 'temperatur

28、e' fit:%X In put : x%Y Output: y% Output:% fitresult : a fit object represe nting the fit.% gof : structure with good ness-of fit info. % 另请参阅 FIT, CFIT, SFIT.% 由 MATLAB 于 17-Aug-2015 19:27:00 自动生成% Fit: 'temperature'.xData, yData = prepareCurveData( x, y );% Set up fittype and options.f

29、t = fittype( 'fourier2');opts = fitopti on s( 'Method', 'Non li nearLeastSquares'); opts.Display = 'Off;opts.Normalize = 'on'opts.StartPoi nt = 0 0 0 0 0 1.9316882339819;% Fit model to data.fitresult, gof = fit( xData, yData, ft, opts );% Plot fit with data.figure

30、( 'Name', 'temperature');h = plot( fitresult, xData, yData );lege nd( h, 'y vs. x', 'temperature', 'Locati on', 'NorthEast'); % Label axesxlabel( 'x');ylabel( 'y');grid onnihe.mx=1:24;y=12 12 13 13 13 13 13 13 12 13 13 15 17 18 19 21 22

31、 22 21 20 16 13 12 12; fitresult, gof = createFit(x, y)偏微分方程求解程序function U x t=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C)%古典显式格式求解抛物型偏微分方程%U x t=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C)%方程：u_t=C*u_xx 0 <= x <= uX,0<= t <= uT% 初值条件：u(x,0)=phi(x)%边值条件:u(0,t)=p

32、si1(t), u(uX,t)=psi2(t)%输出参数：U -解矩阵，第一行表示初值，第一列和最后一列表示边值，第二 行表示第2层%x -空间变量%t -时间变量%输入参数：uX -空间变量x的取值上限%uT -时间变量t的取值上限%phi -初值条件，定义为内联函数%psi1 -边值条件，定义为内联函数%psi2 -边值条件，定义为内联函数%M -沿x轴的等分区间数%N -沿t轴的等分区间数%C -系数，默认情况下C=1%应用举例：% uX=1;uT=0.2;% M=15;N=100;C=1;% phi=i nlin e('si n( pi*x)');psi1=i nli

33、ne('0');psi2=i nlin e('O');%U x t=PDEParabolicClassicalExplicit(uX,uT,phi,psi1,psi2,M,N,C);%inlin e('263+15.11+2.216*cos(1.757*t)+3.687*si n(1.757*t)-1.895*cos(2*1.757*t)+1. 728*si n(2*1.757*t)')%设置参数C的默认值if nargin=7C=1;end %计算步长dx=uX/M;%x 的步长 dt=uT/N;%t 的步长 x=(O:M)*dx;t=(O:N

34、)*dt;r=C*dt/dx/dx;% 步长比r1=1-2*r;if r > 0.5disp('r > 0.5,不稳定')end%计算初值和边值U=zeros(M+1,N+1);for i=1:M+1U(i,1)=phi(x(i);endfor j=1:N+1U(1,j)=psi1(t(j);U(M+1,j)=psi2(t(j);end%逐层求解for j=1:Nfor i=2:MU(i,j+1)=r*U(i-1,j)+r1*U(i,j)+r*U(i+1,j); endendU=U'%作出图形mesh(x,t,U);title('温度分布图像'

35、;)ylabel('空间变量x')xlabel('时间变量t ')zlabel('温度 U')return;18温度随时间，随距离变化的图像19#°0 102 J 30 A 05D.6 D7 0.80.9诃间娈童1#图7第二层温度随时间的变化图像#3420101 H I? 1II- I 1 nn 09on? orw 0 06 0090*012 0.1J 0 5 01002时问直童t图9第二层温度随时间的变化图像空问叢悔图10第三层温度随距离的变化图像21#阖庫分布凶呀河间噜童M#图11第四层温度随时间的变化图像品宸分布图備图12第四层温

36、度随距离的变化图像问题二模型二误差分析程序：%求土壤温度分布clear;clc;format short ea=input('请输入系数a的值:');l=input('请输入长度l的值：)；M=input('请输入将区间0,1等分的个数M :');ot=input('请输入时间增量ot的值：)； n=input('请输入运行次数n的值：)； ox=l/M;x0=zeros(M+1,1);for ii=1:Mx0(ii+1)=ii*ox;endu=sin(pi*x0/l);%t=0 时 u(x,t)的值 r=aA2*ot/(ox)A2;fo

37、r ii=1: n%数据的输入B=zeros(M-1,1);%存放系数矩阵主对角线元素A=zeros (M-2,1);%存放系数矩阵主对角线元素下方次对角线的元素 C=zeros (M-2,1);%存放系数矩阵主对角线元素上方次对角线的元素 S=zeros(M-1,1);%存放右端的常数项for ii=1:M-2B(ii)=1+2*r;A(ii)=-r;C(ii)=-r;S(ii)=u(ii+1,1);endB(M-1)=1+2*r;S(M-1)=u(M,1);u(1,2)=0;u(M+1,2)=0;S(1,1)=S(1,1)+r*u(1,2);S(M-1,1)=S(M-1,1)+r*u(M+

38、1,2);%追赶法S(1)=S(1)/B(1);T=B(1);k=2;while k=MB(k-1)=C(k-1)/T;T=B(k)-A(k-1)*B(k-1);S(k)=(S(k)-A(k-1)*S(k-1)/T;k=k+1;endk=1;while k=M-1S(M-1-k)=S(M-1-k)-B(M-1-k)*S(M-k);k=k+1;endu(2:M,2)=S; %把结果放入矩阵u中u(:,1)=u(:,2);% 过河拆桥end%计算精确值，存放在u的第二列for x=0:Mu(x+1,2)=exp(-(pi*a/l)A2* n*ot)*si n(pi*x*ox/l);end%计算最大

39、相对误差ez=zeros(M-1,1);for ii=2:Mez(ii-1)=abs(u(ii,1)-u(ii,2)/u(ii,2); endE=max(ez);fprintf ('最后时刻数值解与精确解分别为:n');disp(u);fprintf ('差分法得到的结果与正确结果的最大相对误差为：');disp(num2str(E*100)'%');%画二维图比较plot(x0,u(:,1),'g-',x0,u(:,2),'m*');legend('数值解','精确解')xlabel('x'),ylabel('u(x,t)')title('数值解与精确解比较')问题三温度含水量拟合程序creatFit.m function fitresult, gof = createFit(x, y)%CREATEFIT(X,Y)% Create a fit.% Data for 'untitled fit 1' fit:%X In put : x%Y Output: y% Output:% fitresult : a fit object