2016年1月29日到2月2日两个计算机系的和一个数学系的大神开始了他们的愉快的数学建模美赛之旅。今年的MCM有三道题，A题（连续型）是A Hot Bath，B题（离散型）是Space Junk，C题（大数据型）是The Goodgrant Challenge。经过一天的思考和讨论，发现A题比较贴近实际，需要用到热力学和动力学知识；B题需要查找大量的数据，翻阅大量的文献，同时需要用到航天物理知识；C题数据可以直接导入数据库，但是项目类别较多，而且数据挖掘算法需要自己思考，虽然我是计算机专业的，我也不想选这题……最后就敲定做A题了。
第二天看了两眼热力学和流体力学的书，整个人都是不好的……面对满书的公式，感觉自己并不能帮上什么忙。好在数学系的同学很快给我一个公式，说可以用MATLAB进行计算。于是兴冲冲的下了个MATLAB R2014b。
先说说这个三维热传导公式：
于是我在网上搜MATLAB的偏微分方程的相关内容。方法有很多：
- 使用图形化界面的PDE Tool，但是只能接二维特殊形式的偏微分方程；
- 使用pdepe()函数，可以求解一微一般形式的偏微分方程；
- 使用elliptic()，parabolic()，hyperbolic()函数求解特殊形式的偏微分方程。
经过确认，要求解的偏微分方程为抛物线型。又在网上找了半天，都说第三种方法只能解二维的偏微分方程。就在准备要放弃的时候，我想到了可以查官方文档。官方是这么写的：
那就说明可以解三维的！
先试试官方的例子，但是第一句话都没法正常运行。
model = createpde();
不应该啊，这可是官方文档。排查了半天，可能是MATLAB版本的问题。官方文档应该是提供给最新的R2015b的，这也许是R2015b新加入的函数呢？于是抓紧时间下载并安装了最新的R2015b（这其中还有一段小插曲），然后运行官方的例子，就没有问题了。下面附上代码及简单说明。
%新建模型（一个立方体中间刨去一个圆柱） [xg, yg] = meshgrid(-3:0.25:3,-1.5:0.25:1.5); xg = xg(:); yg = yg(:); t = (pi/24:pi/24:2*pi)'; x = cos(t); y = sin(t); circShp = alphaShape(x,y,2); in = inShape(circShp,xg,yg); xg = [xg(~in); cos(t)]; yg = [yg(~in); sin(t)]; zg = ones(numel(xg),1); xg = repmat(xg,5,1); yg = repmat(yg,5,1); zg = zg*(0:.25:1); zg = zg(:); shp = alphaShape(xg,yg,zg); [elements,nodes] = boundaryFacets(shp); nodes = nodes'; elements = elements'; model = createpde(); geometryFromMesh(model,nodes,elements); generateMesh(model); %显示所建模型的图形 handl = pdegplot(model,'FaceLabels','on'); view(-42,24) handl(1).FaceAlpha = 0.5; %添加边界条件，设置1号面恒为37度（狄利克雷条件） applyBoundaryCondition(model,'Face',1,'u',37); generateMesh(model); p = model.Mesh.Nodes; x = p(1,:); y = p(2,:); z = p(3,:); %设置初始温度条件 u0 = 42-0.5*(x-3).*(x-3)-0.5*(y-1.5).*(y-1.5)-0.5*(z-1).*(z-1); %设置参数 d = 4.2*10^6; c = 5000; a = 1000; f = 1000*25; [Kc,Fc,B,ud] = assempde(model,c,a,f); [~,M,~] = assema(model,0,d,f); %设置计算区间，t从0到1800秒，每60秒算一次结果 tlist = linspace(0,1800,31); u = parabolic(u0,tlist,Kc,Fc,B,ud,M); %计算最小温度和最大温度 umin = min(min(u)); umax = max(max(u)); %显示第一个子图，显示第1个结果 subplot(2,2,1) pdeplot3D(model,'colormapdata',u(:,1)) colorbar view(125,22) title 't = 0' caxis([umin umax]); %显示第二个子图，显示第11个结果 subplot(2,2,2) pdeplot3D(model,'colormapdata',u(:,11)) colorbar view(125,22) title 't = 10' caxis([umin umax]); %显示第三个子图，显示第21个结果 subplot(2,2,3) pdeplot3D(model,'colormapdata',u(:,21)) colorbar view(125,22) title 't = 20' caxis([umin umax]); %显示第四个子图，显示第31个结果 subplot(2,2,4) pdeplot3D(model,'colormapdata',u(:,31)) colorbar view(125,22) title 't = 30' caxis([umin umax]);
运行的结果如图所示：
附：比赛最后一天真的是越写越high……然后五点交完论文滚去睡觉了。睡到十二点起来聚餐。这次美赛真的是很有意思（还是队友比较有意思→_→）。