一些关于CUMCM的思考

去年参加的CUMCM，觉得代码闲置在附录里边比较可惜就打算分享出来了。
国赛的题目想要求出结果并不是很难，使得论文的写作以及验证等环节愈发重要，语言表达是否精练以及模型的评价与改进是否合理有效成为了评价参赛作品的一项重要指标。由于去年作品要求通过客户端上传，一定要注意上传支撑材料的大小不要超过20M，注意参赛时间安排以及qq账号的清洁（！）。
本文给出的思路，旨在辅助对于题目的理解并帮助解释贴出的Matlab代码，不代表最优做法。本文思路及代码仅供参考，思路中涉及到的公式、方法以及详细步骤等请移步至参考文献。代码中变量命名比较乱，仅供参考。

CUMCM2018 A题

A题 高温作业专用服装设计
在高温环境下工作时，人们需要穿着专用服装以避免灼伤。专用服装通常由三层织物材料构成，记为I、II、III层，其中I层与外界环境接触，III层与皮肤之间还存在空隙，将此空隙记为IV层。
为设计专用服装，将体内温度控制在37ºC的假人放置在实验室的高温环境中，测量假人皮肤外侧的温度。为了降低研发成本、缩短研发周期，请你们利用数学模型来确定假人皮肤外侧的温度变化情况，并解决以下问题：
(1) 专用服装材料的某些参数值由附件1给出，对环境温度为75ºC、II层厚度为6 mm、IV层厚度为5 mm、工作时间为90分钟的情形开展实验，测量得到假人皮肤外侧的温度（见附件2）。建立数学模型，计算温度分布，并生成温度分布的Excel文件（文件名为problem1.xlsx）。
(2) 当环境温度为65ºC、IV层的厚度为5.5 mm时，确定II层的最优厚度，确保工作60分钟时，假人皮肤外侧温度不超过47ºC，且超过44ºC的时间不超过5分钟。
(3) 当环境温度为80 时，确定II层和IV层的最优厚度，确保工作30分钟时，假人皮肤外侧温度不超过47ºC，且超过44ºC的时间不超过5分钟。

附件1. 专用服装材料的参数值
附件2. 假人皮肤外侧的测量温度

对于题目的分析

对于该题目可以建立如下假设：
1.假设外界环境与假人体内的气温保持不变。
2.假设辐射传热的影响可以忽略。
3.假设人的体型对传热的影响可以忽略，人体可以等效为圆柱体。
4.假设热量在材料交接交界面处的传导方向为交界面法线方向。

这是那种第一问做出来之后第二问第三问都迎刃而解的题目，找准建立模型的物理理论依据，再确保了解决偏微分方程的求解手段就可以得到答案了。

第一问

1.1 总体思路

问题1要求在给定外界温度和高温工作服各层材料厚度的条件下求解工作服的温度分布情况。根据热学分析可以建立出关于各层材料传热情况的偏微分方程组模型。工作服的初始温度以及工作过程中材料各交界面的温度作为各层偏微分方程边界条件。由于此微分方程组模型较为复杂，运用有限差分法将微分方程化为差分方程，再将差分方程组以矩阵的形式写出，并用追赶法（Thomas算法）求解差分方程组模型可得出工作服的温度分布。
针对微分的边界条件确定，最内层材料与假人皮肤交界面的温度同最外层材料与外界空气交界面的温度依据牛顿冷却定律通过外界环境气温、假人内部温度计算，运用牛顿冷却定律时需要确定对流换热系数这一未知参数，此外工作服的初始温度也是未知量。对于这两个未知量，先大致估计这两个未知量的范围，再按照一定步长遍历搜索。对于每个遍历到的值，代入模型中运用有限差分法计算，得出相应的皮肤表面的温度，与附件二给出的实测值进行对比，取残差和最小时的值作为初始温度和对流换热系数的理想值。

1.2 模型求解

1.2.1 偏微分方程的建立
热传导在不同温度的物体相互接触或者同一物体不同部分具有不同温度时产生，符合本题中几种材料互相接触传递热量的情形。将假人抽象为多层互相包裹的同轴圆柱体后，可以对上述三维几何模型进行降维处理，建立一维几何模型。且该题可忽略对流传热与辐射传热的影响，物体内部的热传递由傅里叶定律确定，而物体表面和内部的热量流动又满足能量守恒定律，依据题目给定部分边界条件可建立一维热传导模型。
物体在边界上与其他传热介质接触时，物体表面与介质的温度往往不相同，此时运用牛顿冷却定律：单位时间从物体流入介质中的热量和两者的温度差成正比，可确定热传导的第三类边界条件。
1.2.2 有限差分法——化为差分方程
生成网格–>微分化差分–>差分方程矩阵化–>边界条件矩阵化–>合成矩阵差分方程组（这一步比较考验耐心）–>追赶法
具体做法可参考这篇文章：史策.热传导方程有限差分法的MATLAB实现[J].咸阳师范学院学报,2009,24(04):27-29+36.
1.2.3 追赶法——解差分方程
实际建立出方程组的系数矩阵为三对角矩阵，即非零元素都集中在主对角线和其相邻两条次对角线上，故容易分解成两个三角阵，因此可以用追赶法求解。
追赶法的理论可参考这篇博文：追赶法求解三对角方程组
1.2.4 输出结果至Excel文件。

1.3 Matlab实现

代码如下：

%%为应用牛顿冷却定律确定边值条件，需要在测试范围内遍历k21取值，计算边值条件与附件二中数据进行比对；
%%选取得到最小平均残差，并将最大残差控制在一定范围内的k21值（此程序得到k21=-0.0046）；
clc;
clear;
load('tX4.mat');%导入附件2假人皮肤外侧温度变化数据，事先已存储为变量，在此问作为参考值使用
%%常量及边界条件定义
d1=0.6; d2=6; d3=3.6; d4=5;%单位为mm
k1=0.082; k2=0.370; k3=0.045; k4=0.028;%热传导率
c1=1377; c2=2100; c3=1726; c4=1005;%比热
p1=300; p2=862; p3=74.2; p4=1.18;%密度
a2=[k1/(c1*p1) k2/(c2*p2) k3/(c3*p3) k4/(c4*p4)];%a2=k/(cp)
dx=0.1;
m=floor([d1/dx d2/dx d3/dx d4/dx]);
m1=m(1);m2=m(2);m3=m(3);m4=m(4);
for k=1:(m1+m2+m3+m4)u(k,1)=14.7;%设衣服初始温度为14.7度(最小平均残差)、25度(室温)、37度(假人皮肤外温度)时分别进行残差计算
end
N=5400;%离散化点数k11=-1.2;k21=-0.0046;
for j=1:NM=m1+m2+m3+m4;ox1=d1/m1;ox2=d2/m2;ox3=d3/m3;ox4=d4/m4;ot=1;%%变量定义r1=a2(1)*ot/(ox1/1000)^2;r2=a2(2)*ot/(ox2/1000)^2;r3=a2(3)*ot/(ox3/1000)^2;r4=a2(4)*ot/(ox4/1000)^2;%数据的输入B=zeros(M-1,1);%存放系数矩阵主对角线元素A=zeros(M-1,1);%存放系数矩阵主对角线元素下方次对角线的元素C=zeros(M-1,1);%存放系数矩阵主对角线元素上方次对角线的元素S=zeros(M-1,1);%存放右端的常数项for ii=1:m1B(ii)=1+2*r1;A(ii)=-r1;C(ii)=-r1;S(ii)=u(ii+1,1);endfor ii=m1+1:m1+m2B(ii)=1+2*r2;A(ii)=-r2;C(ii)=-r2;S(ii)=u(ii+1,1);endfor ii=m1+m2+1:m1+m2+m3B(ii)=1+2*r3;A(ii)=-r3;C(ii)=-r3;S(ii)=u(ii+1,1);endfor ii=m1+m2+m3+1:m1+m2+m3+m4-1B(ii)=1+2*r4;A(ii)=-r4;C(ii)=-r4;S(ii)=u(ii+1,1);endA(1)=0;C(M-1)=0;B(M-1)=1+2*r4;S(M-1)=u(M,1);datac(j)=(u(2,1)-dx*75*k11/k1)/(1-dx*k11/k1);datad(j)=(dx*37*k21/k4-u(end-1,1))/(dx*k21/k4-1);%根据牛顿冷却定律计算出人皮肤外侧温度作为边界条件u(1,2)=datac(j);%热源边值条件输入！！u(M+1,2)=datad(j);%人皮肤表面温度：边值条件输入！！S(1,1)=r1*u(1,2)+S(1,1);S(M-1,1)=S(M-1,1)+r4*u(M+1,2);%…矩阵建造结束%追赶法求解矩阵方程，b是对角线元素，a为-1对角线，c为1对角线T=B(1);U=B;l=B;for k=2:M-1  %解Lx=b：追l(k)=A(k)/U(k-1);U(k)=B(k)-C(k-1)*l(k);endY=S;for k=2:M-1Y(k)=S(k)-l(k)*Y(k-1);%中间解endX=Y;X(M-1)=Y(M-1)/U(M-1);for k=M-2:-1:1  %解Ux=b：赶X(k)=(Y(k)-C(k)*X(k+1))/U(k);endoutputa(j,:)=X;u(2:M,2)=X; %把结果放入矩阵 u 中u(:,1)=u(:,2);% 计算下一时刻
end
dada=max(abs(outputa(31:end,end)-tX4(32:end,2)))%最大残差
kaka=mean(abs(outputa(31:end,end)-tX4(32:end,2)))%平均残差

clc;
clear;
load('tX4.mat');%导入附件2假人皮肤外侧温度变化数据，事先已存储为变量
%%常量定义
d1=0.6; d2=6; d3=3.6; d4=5;%各段厚度，单位为mm
k1=0.082; k2=0.370; k3=0.045; k4=0.028;%热传导率
c1=1377; c2=2100; c3=1726; c4=1005;%比热
p1=300; p2=862; p3=74.2; p4=1.18;%密度
a2=[k1/(c1*p1) k2/(c2*p2) k3/(c3*p3) k4/(c4*p4)];%a2=k/(cp)
dx=0.01;%离散化步长
m=floor([d1/dx d2/dx d3/dx d4/dx]);
m1=m(1);m2=m(2);m3=m(3);m4=m(4);
for k=1:(m1+m2+m3+m4)u(k,1)=37;
end
N=5400;%离散化点数
for j=1:NM=m1+m2+m3+m4;ox1=d1/m1;ox2=d2/m2;ox3=d3/m3;ox4=d4/m4;ot=1;%%变量定义r1=a2(1)*ot/(ox1/1000)^2;r2=a2(2)*ot/(ox2/1000)^2;r3=a2(3)*ot/(ox3/1000)^2;r4=a2(4)*ot/(ox4/1000)^2;%数据的输入B=zeros(M-1,1);%存放系数矩阵主对角线元素A=zeros(M-1,1);%存放系数矩阵主对角线元素下方次对角线的元素C=zeros(M-1,1);%存放系数矩阵主对角线元素上方次对角线的元素S=zeros(M-1,1);%存放右端的常数项for ii=1:m1B(ii)=1+2*r1;A(ii)=-r1;C(ii)=-r1;S(ii)=u(ii+1,1);endfor ii=m1+1:m1+m2B(ii)=1+2*r2;A(ii)=-r2;C(ii)=-r2;S(ii)=u(ii+1,1);endfor ii=m1+m2+1:m1+m2+m3B(ii)=1+2*r3;A(ii)=-r3;C(ii)=-r3;S(ii)=u(ii+1,1);endfor ii=m1+m2+m3+1:m1+m2+m3+m4-1B(ii)=1+2*r4;A(ii)=-r4;C(ii)=-r4;S(ii)=u(ii+1,1);endk11=-0.9;A(1)=0;C(M-1)=0;B(M-1)=1+2*r4;S(M-1)=u(M,1);datac(j)=(u(2,1)-dx*75*k11/k1)/(1-dx*k11/k1);u(1,2)=datac(j);%热源边界条件输入！！u(M+1,2)=tX4(j,2);%人皮肤表面温度：边值条件输入！！S(1,1)=r1*u(1,2)+S(1,1);S(M-1,1)=S(M-1,1)+r4*u(M+1,2);%…矩阵建造结束%追赶法求解矩阵方程，b是对角线元素，a为-1对角线，c为1对角线T=B(1);U=B;l=B;for k=2:M-1  %解Lx=b：追l(k)=A(k)/U(k-1);U(k)=B(k)-C(k-1)*l(k);endY=S;for k=2:M-1Y(k)=S(k)-l(k)*Y(k-1);%中间解endX=Y;X(M-1)=Y(M-1)/U(M-1);for k=M-2:-1:1  %解Ux=b：赶X(k)=(Y(k)-C(k)*X(k+1))/U(k);endoutputa(j,:)=X;u(2:M,2)=X; %把结果放入矩阵 u 中u(:,1)=u(:,2);% 用来计算下一时刻
end
xlswrite('problem1.xlsx',outputa);%生成Excel文件存储温度分布矩阵，其中第(i,j)个单元格的值代表第i点厚度、第j点时刻的温度值

第二问

2.1 总体思路

问题2是在给定了环境温度和除Ⅱ层之外每一层的厚度的条件下，要求确定Ⅱ层的最优厚度，使得该高温工作服能够满足在65℃环境中工作60分钟时假人皮肤外侧温度不超过47℃且超过44℃的时间不超过5分钟。已知材料越厚则隔绝高温的效果越好，但从实际出发，高温工作服越薄则质量越轻，活动越灵活，就越有利于提高使用者穿着舒适度和正常工作能力。同时，相同材质的高温工作服越薄则制造成本也就越低。本题就是要求出在高温工作服能够达到给定的高温防护标准的前提下Ⅱ层的最小厚度。
可以类比问题1的做法建立列写热传导微分方程建立模型。本问中皮肤外侧的温度未知，只知道假人体内温度恒定为37℃，因此在确定假人皮肤外侧与Ⅳ层交界面的温度这一边界条件时需要用到牛顿冷却定律，这与问题1中确定最外层材料与空气交界面温度以及假人皮肤外侧与Ⅳ层交界面温度的做法一致。

2.2 模型求解

对于建立好的微分方程组模型，沿用问题1的有限差分法进行求解。为了找到满足要求的最小Ⅱ层厚度，可先选定一个较小的Ⅱ层厚度的值，计算假人皮肤外侧的温度随时间的变化情况，遍历搜索的算法逐渐增大Ⅱ层材料厚度的值，直至皮肤外侧的温度刚好满足题目要求，此时的Ⅱ层材料厚度即为最佳厚度。

2.3 Matlab实现

代码如下：

%%最优厚度即满足要求的最小厚度，此程序选择满足问题二要求的第二层最小厚度；
clc;
clear;
%%常量及边界条件定义
for dd=8.5:0.01:9.5  %遍历范围1：5~10，步长：0.1；遍历范围2:8.5~9.5，步长：0.01；
d1=0.6; d2=dd; d3=3.6; d4=5.5;%单位为mm
k1=0.082; k2=0.370; k3=0.045; k4=0.028;%热传导率
c1=1377; c2=2100; c3=1726; c4=1005;%比热
p1=300; p2=862; p3=74.2; p4=1.18;%密度
a2=[k1/(c1*p1) k2/(c2*p2) k3/(c3*p3) k4/(c4*p4)];%a2=k/(cp)
dx=0.1;
m=floor([d1/dx d2/dx d3/dx d4/dx]);
m1=m(1);m2=m(2);m3=m(3);m4=m(4);
for k=1:(m1+m2+m3+m4)u(k,1)=14.7;%由于此温度时模拟边值条件平均残差最小，在第二问遍历寻优时利用该温度作为服装的初始温度
end
N=3600;%离散化点数
outputa=zeros(3600,m1+m2+m3+m4-1);k11=-1.2;
for j=1:NM=m1+m2+m3+m4;ox1=d1/m1;ox2=d2/m2;ox3=d3/m3;ox4=d4/m4;ot=1;%%变量定义r1=a2(1)*ot/(ox1/1000)^2;r2=a2(2)*ot/(ox2/1000)^2;r3=a2(3)*ot/(ox3/1000)^2;r4=a2(4)*ot/(ox4/1000)^2;%数据的输入B=zeros(M-1,1);%存放系数矩阵主对角线元素A=zeros(M-1,1);%存放系数矩阵主对角线元素下方次对角线的元素C=zeros(M-1,1);%存放系数矩阵主对角线元素上方次对角线的元素S=zeros(M-1,1);%存放右端的常数项for ii=1:m1B(ii)=1+2*r1;A(ii)=-r1;C(ii)=-r1;S(ii)=u(ii+1,1);endfor ii=m1+1:m1+m2B(ii)=1+2*r2;A(ii)=-r2;C(ii)=-r2;S(ii)=u(ii+1,1);endfor ii=m1+m2+1:m1+m2+m3B(ii)=1+2*r3;A(ii)=-r3;C(ii)=-r3;S(ii)=u(ii+1,1);endfor ii=m1+m2+m3+1:m1+m2+m3+m4-1B(ii)=1+2*r4;A(ii)=-r4;C(ii)=-r4;S(ii)=u(ii+1,1);endk21=-0.0046;A(1)=0;C(M-1)=0;B(M-1)=1+2*r4;S(M-1)=u(M,1);datac(j)=(u(2,1)-dx*65*k11/k1)/(1-dx*k11/k1);datad(j)=(dx*37*k21/k4-u(end-1,1))/(dx*k21/k4-1);u(1,2)=datac(j);%热源边值条件输入！！u(M+1,2)=datad(j);%人皮肤表面温度：边值条件输入！！S(1,1)=r1*u(1,2)+S(1,1);S(M-1,1)=S(M-1,1)+r4*u(M+1,2);%…矩阵建造结束%追赶法求解矩阵方程，b是对角线元素，a下c上T=B(1);U=B;l=B;for k=2:M-1  %解Lx=b：追sl(k)=A(k)/U(k-1);U(k)=B(k)-C(k-1)*l(k);endY=S;for k=2:M-1Y(k)=S(k)-l(k)*Y(k-1);%中间解endX=Y;X(M-1)=Y(M-1)/U(M-1);for k=M-2:-1:1  %解Ux=b：赶X(k)=(Y(k)-C(k)*X(k+1))/U(k);endoutputa(j,:)=X;u(2:M,2)=X; %把结果放入矩阵 u 中u(:,1)=u(:,2);% 下一时刻的值
end
if outputa(3300-1,end)>=44;%条件1continue;
end
if outputa(:,end)>=47;%条件2continue;
end
outi=dd;
break;%得到最小厚度即退出循环
end

第三问

3.1 总体思路

问题3与问题2类似，也是在给定环境气温和隔热效果要求的条件下求高温工作服某些层的厚度。与问题2不同的是，本问要同时确定Ⅱ、Ⅳ两层材料的最优厚度。此处的厚度最优一方面要使得两种材料都尽可能薄，另一方面，考虑到第Ⅱ层材料为某种制作高温工作服的原材料，第Ⅳ层为衣服与皮肤之间的空气间隙，从附件给出的参数上看，空气的密度远远小于Ⅱ层材料的密度，因此同等厚度的Ⅱ层材料的质量远远大于Ⅳ层空气的质量，而且Ⅱ层材料需要成本，而Ⅳ层空气间隙不需成本，所以无论从工作服舒适度还是成本上看，在进行设计时，应该优先使Ⅱ层材料的厚度尽量小，其次再考虑减小Ⅳ层空气间隙的厚度。
在问题2的求解过程中已经求得了Ⅰ层材料与外界空气的热交换系数以及Ⅳ层材料与假人皮肤的热交换系数，并且已经查阅文献得知在材料与介质不变时热交换系数可认为是定值，因此问题二中求出的值在问题3中可以直接使用。
反复改变Ⅱ、Ⅳ两层材料的厚度值，在使得假人皮肤表面温度满足题目要求的解中寻找出最优的Ⅱ、Ⅳ两层材料的厚度值。
问题3与问题2相比，热量传递的物理模型相同，仅仅有参数上的差别，因此可以直接沿用问题2中的热传导方程建立偏微分方程组模型。

3.2 模型求解

由于本问题的微分方程组模型与问题2形式一致，因此问题3也沿用问题1和问题2求解微分方程组的办法，即用有限差分法求解此模型。在前两问中已经求得热交换系数的值，考虑到热交换系数的值在材料以及介质相同时变化较小，因此本问在求解时直接沿用问题2中的热传导系数值。
为了得到第Ⅱ、Ⅳ两层的最优厚度，使用在一定范围内遍历搜索的算法。最优厚度是要使得两层材料都尽可能薄，并且优先保证第Ⅱ层材料的厚度最小。做两层循环就好。

3.3 Matlab实现

代码如下：

clc;
clear;
%%常量及边界条件定义
num=0;
for dd=8:0.1:15flag=0;for dd1=1.4:0.1:6.4
d1=0.6; d2=dd; d3=3.6; d4=dd1;%单位为mm
k1=0.082; k2=0.370; k3=0.045; k4=0.028;%热传导率
c1=1377; c2=2100; c3=1726; c4=1005;%比热
p1=300; p2=862; p3=74.2; p4=1.18;%密度
a2=[k1/(c1*p1) k2/(c2*p2) k3/(c3*p3) k4/(c4*p4)];%a2=k/(cp)
dx=0.1;
m=floor([d1/dx d2/dx d3/dx d4/dx]);
m1=m(1);m2=m(2);m3=m(3);m4=m(4);
for k=1:(m1+m2+m3+m4)u(k,1)=14.7;
end
N=1800;
outputa=zeros(1800,m1+m2+m3+m4-1);k11=-1.2;
for j=1:NM=m1+m2+m3+m4;ox1=d1/m1;ox2=d2/m2;ox3=d3/m3;ox4=d4/m4;ot=1;%%变量定义r1=a2(1)*ot/(ox1/1000)^2;r2=a2(2)*ot/(ox2/1000)^2;r3=a2(3)*ot/(ox3/1000)^2;r4=a2(4)*ot/(ox4/1000)^2;%数据的输入B=zeros(M-1,1);%存放系数矩阵主对角线元素A=zeros(M-1,1);%存放系数矩阵主对角线元素下方次对角线的元素C=zeros(M-1,1);%存放系数矩阵主对角线元素上方次对角线的元素S=zeros(M-1,1);%存放右端的常数项for ii=1:m1B(ii)=1+2*r1;A(ii)=-r1;C(ii)=-r1;S(ii)=u(ii+1,1);endfor ii=m1+1:m1+m2B(ii)=1+2*r2;A(ii)=-r2;C(ii)=-r2;S(ii)=u(ii+1,1);endfor ii=m1+m2+1:m1+m2+m3B(ii)=1+2*r3;A(ii)=-r3;C(ii)=-r3;S(ii)=u(ii+1,1);endfor ii=m1+m2+m3+1:m1+m2+m3+m4-1B(ii)=1+2*r4;A(ii)=-r4;C(ii)=-r4;S(ii)=u(ii+1,1);endk21=-0.0046;A(1)=0;C(M-1)=0;B(M-1)=1+2*r4;S(M-1)=u(M,1);datac(j)=(u(2,1)-dx*80*k11/k1)/(1-dx*k11/k1);datad(j)=(dx*37*k21/k4-u(end-1,1))/(dx*k21/k4-1);u(1,2)=datac(j);%热源边值条件输入！！u(M+1,2)=datad(j);%人皮肤表面温度：边值条件输入！！S(1,1)=r1*u(1,2)+S(1,1);S(M-1,1)=S(M-1,1)+r4*u(M+1,2);%…矩阵建造结束%追赶法求解矩阵方程，b是对角线元素，a为-1对角线，c为1对角线T=B(1);U=B;l=B;for k=2:M-1  %解Lx=b：追l(k)=A(k)/U(k-1);U(k)=B(k)-C(k-1)*l(k);endY=S;for k=2:M-1Y(k)=S(k)-l(k)*Y(k-1);%中间解endX=Y;X(M-1)=Y(M-1)/U(M-1);for k=M-2:-1:1  %解Ux=b：赶X(k)=(Y(k)-C(k)*X(k+1))/U(k);endoutputa(j,:)=X;u(2:M,2)=X; %把结果放入矩阵 u 中u(:,1)=u(:,2);% 下一时刻值
end
if outputa(1500-1,end)>=44  %条件1continue;
end
if outputa(:,end)>=47  %条件2continue;
end
num=num+1;
outi(num)=dd;
outj(num)=dd1;
if (outputa(1500-1,end)<44)&(outputa(:,end)<47)flag=1;break;%得到每个满足要求的第二层厚度对应的第四层最小厚度即退出循环
endend
if flag==1break;
end
end

另一种做法

还有一种不应用牛顿冷却定律确定人体皮肤外表面温度这个边界条件，而假想出“第五层材料”的做法，虽然可以和本文提供的方法得到一样的结果（第二问、第三问），但是略微欠缺一些物理理论的支持。
如果只是想看到第一问的分布图的话，数值仿真软件Comsol的热传递模块可以较为轻松地得到第一问的仿真结果，但是内存消耗比较严重不是很友好。
——————————————————————————————————————
本文给出的思路，旨在辅助对于题目的理解并帮助解释贴出的Matlab代码，不代表最优做法。本文思路及代码仅供参考，思路中涉及到的公式、方法以及详细步骤等请移步至参考文献。代码中变量命名比较乱，仅供参考。

参考文献

[1]王元明.数学物理方程与特殊函数[M].高等教育出版社.2011,:84-97.
[2] (美) 施罗德 (Schroeder,D.V.) , 著. 热物理学导论[M].北京:世界图书出版公司北京公司, 2008.
[3]史策.热传导方程有限差分法的MATLAB实现[J].咸阳师范学院学报,2009,24(04):27-29+36.
[4]南京大学数学系计算数学专业.偏微分方程数值解法[M].北京.科学出版社.1979