matlab的v4插值
griddata(x,y,z,linspace(min(x),max(x))’,linspace(min(y),max(y))


根据帮助文件，v4插值为双调和样条插值（biharmonic spline interpolation）




关于双调和样条插值，参考的此网址中的文献

http://www.fon.hum.uva.nl/praat/manual/biharmonic_spline_interpolation.html


文献（doi：10.1007/s11004-011-9346-5）中给出的方法：




根据文献所给方法得到的matlab代码（参数仅适用于此案例，具体情况应按文献修改）


   x=[27.5	77.5	111.25	0	0	0	-27.5	-77.5	-111.25	0	0	0];
    y=[0	0	0	27.5	77.5	111.25	0	0	0	-27.5	-77.5	-111.25];
    z=[0.256700000000000,0.0720900000000000,0.0357500000000000,0.246100000000000,0.273500000000000,0.0820800000000000,0.257700000000000,0.0710900000000000,0.0357500000000000,0.236100000000000,0.283500000000000,0.0820800000000000];
   %[X,Y,Z]=griddata(x,y,z,linspace(min(x),max(x))',linspace(min(y),max(y)),'v4');%插值
%Z(find(Z<0))=0;
%以下内容为根据文献得出代码，作用相当于上方两行代码
    X=linspace(min(x),max(x));%将最小值与最大值划分为100块
    Y=linspace(min(y),max(y));
 %以下求权重矩阵   
i=1;
while i~=13
    j=1;
    while j~=13
        r(i,j)=((x(i)-x(j))^2+(y(i)-y(j))^2)^0.5;
        if (r(i,j)==0)
            d(i,j)=0;
        else
            d(i,j)=(log(r(i,j))-1)*r(i,j)^2;%格林函数矩阵，即文献中d(i,j)，即G
        end
        j=j+1;
    end
    i=i+1;
end
w=inv(d)*z';%权重矩阵

%以下求XY的100*100矩阵
n=1;
while n~=101
    m=1;
    while m~=101
        X(m,n)=X(1,n);
        
        m=m+1;
    end
    n=n+1;
end

Y=X';

%以下求Z矩阵
n=1;
while n~=101
    m=1;
    while m~=101
        i=1;
        while i~=13
            v(1,i)=((X(m,n)-x(i))^2+(Y(m,n)-y(i))^2)^0.5;%即文献中r(i,j)
            if v(1,i)==0
                u(1,i)=0;
            else
                u(1,i)=(log(v(1,i))-1)*v(1,i)^2;%即文献中d(i,j)，即Gp
            end
            i=i+1;
        end
        Z(m,n)=u*w;
        m=m+1;
    end
    n=n+1;
end
Z(find(Z<0))=0;%将Z中小于0的数取0
figure,surf(X,Y,Z); 



v4插值代码（此结果与标题4的计算结果相同）

x=[27.5	77.5	111.25	0	0	0	-27.5	-77.5	-111.25	0	0	0];
    y=[0	0	0	27.5	77.5	111.25	0	0	0	-27.5	-77.5	-111.25];
    z=[0.256700000000000,0.0720900000000000,0.0357500000000000,0.246100000000000,0.273500000000000,0.0820800000000000,0.257700000000000,0.0710900000000000,0.0357500000000000,0.236100000000000,0.283500000000000,0.0820800000000000];
   [X,Y,Z]=griddata(x,y,z,linspace(min(x),max(x))',linspace(min(y),max(y)),'v4');%插值
Z(find(Z<0))=0;
figure,surf(X,Y,Z); 

6.本人matlab初学者，欢迎批评指正

关于连续信号的傅里叶分析，可以参考博文《学习笔记之——基于matlab的数字通信系统（1）&连续信号的傅里叶分析》

 

目录

离散时间信号的傅里叶变换（DTFT）

连续时间信号的抽样——抽样定理

离散傅里叶级数（DFS）——离散的周期序列的傅里叶级数

离散傅里叶变换（DFT）

快速傅里叶变换FFT

代码

参考文献

 

离散时间信号的傅里叶变换（DTFT）

先给出课本的解析







对于DTFT有如下特点：

1、是离散的，所以变换需要求和而不是积分。

2、是w的连续函数；且是w的周期函数，周期为2pi



3、可以看作将在频域展开为傅里叶级数，傅里叶

DTFT的重要性质：



 

连续时间信号的抽样——抽样定理

实际种，信号的抽样是通过A/D芯片来实现将连续信号变成数字信号。通过抽样，将连续信号的x(t)变为数字信号x(nTs)，将x(t)的傅里叶变换变为





抽样函数为，其对应的傅里叶变换为

对于信号采样，有公式：

根据时域卷积定理，有：

经运算后有：

将连续信号抽样变为离散信号，对应的离散信号的频谱将变成周期的。而变成周期的方法就是进行频谱的周期延拓。如下图所示：



进而得出抽样定理：



离散傅里叶级数（DFS）——离散的周期序列的傅里叶级数



 

离散傅里叶变换（DFT）



在计算机上实现频谱分析及其他方面的处理工作时，对信号的要求是：在时域和频域上都应该是离散的。在四种傅里叶变换种，只有DFS在时域和频域上都说离散的，但是却都是无限长的

对于DFS有



在DFS种，n和k仍取无穷长，但实际上并没有必要，因为他们都是以N点（采样点）为周期的。故此有：



显然，DFT并不是一种新的傅里叶变换，它实际上是来自于DFS，只不过它仅仅在时域和频域各取一个周期而已。

 

快速傅里叶变换FFT

 

代码

DTFT

N=8;                         %原离散信号有8点
n=[0:1:N-1]                  %原信号是1行8列的矩阵
syms t;
xn=0.5.^n;f=0.5.^t                   %构建原始信号，为指数信号
subplot(311);
ezplot(f);
title('原始信号(指数信号)');

subplot(312)
stem(n,xn);
title('原始信号(离散化),');

w=[-800:1:800]*4*pi/800;     %频域共-800----+800 的长度（本应是无穷，高频分量很少，故省去）    
X=xn*exp(-j*(n'*w));         %求dtft变换，采用原始定义的方法，对复指数分量求和而得

subplot(313);
plot(w/pi,abs(X));
title('DTFT变换')



clear all
%%%%%%%%%%%%%%%%
n1=-15:15;
h1=exp(-abs(0.1*n1));
subplot(2,2,1);stem(n1,h1);title('h1');
%%%%%%%%%%%%%%%%%%%%%%%
n2=0:20;
h2(n2+1)=1;
subplot(2,2,2);stem(n2,h2);title('h2');

w=-4:0.001:4;%产生要计算的w的范围
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Hjw1=h1*(exp(-j*pi).^(n1'*w));
subplot(2,2,3);plot(w,abs(Hjw1))
title('H1');xlabel('pi弧度(w)');ylabel('振幅')
%%%%%%%%%%%%%%%%%%%%%%%%%%%
Hjw2=h2*(exp(-j*pi).^(n2'*w));
subplot(2,2,4);plot(w,abs(Hjw2));
title('H2');xlabel('pi弧度(w)');ylabel('振幅')

%%%%%%%非周期离散信号的DTFT是——周期性，且连续的



验证卷积定理



clear all
n=0:30;
h=sinc(0.2*n);
subplot(3,3,1);stem(n,h);title('系统的时域响应')
x=2*sin(0.2*pi*n)+3*cos(0.4*pi*n);
subplot(3,3,2);stem(n,x);title('系统的输入')
%%%%%%%%%%%%%%%%%%%DTFT
w=-1:0.001:1;%%%要计算的频率的范围
Hjw=h*(exp(-j*pi).^(n'*w));
subplot(3,3,4);plot(w,abs(Hjw))
title('系统的频域响应');xlabel('pi弧度(w)');ylabel('振幅')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Xjw=x*(exp(-j*pi).^(n'*w));
subplot(3,3,5);plot(w,abs(Xjw));
title('系统的输入的频域');xlabel('pi弧度(w)');ylabel('振幅')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Yjw=Xjw.*Hjw;
%%%%%%%%%%%%%%%%%
subplot(3,3,6);plot(w,abs(Yjw));
title('系统的输出的频域');xlabel('pi弧度(w)');ylabel('振幅')


n1=0:2*length(n)-2;%%%%%%%%计算IDTFT
dw=0.001*pi;%%%%%%%%%%分段求和来实现积分
y=(dw*Yjw*(exp(j*pi).^(w'*n1)))/(2*pi);
subplot(3,3,7);stem(abs(y));title('通过IDTFT计算出的输出序列Y');

y1=conv(x,h);
subplot(3,3,9);stem(abs(y1));title('通过时域卷积计算出的输出序列Y1')







DFT变换

clear all
n=0:30;%%%%%%%%n的范围
x=sin(0.2*n).*exp(-0.1*n);
subplot(2,1,1);stem(n,x);title('序列x')
k=0:30;
N=31;
Wnk=exp(-j*2*pi/N).^(n'*k);%%%%%%DFT定义
X=x*Wnk;
subplot(2,1,2);stem(-15:15,[abs(X(17:end)) abs(X(1:16))])%%%%%%%%对坐标进行了对称排列
title('X的DFT')




上面的都是直接采用原始的定义公式。在matlab中有FFT函数。

DFT的循环卷积



clear all
%%%%%%产生两个信号序列
h=[6 3 4 2 1 -2];
x=[3 2 6 7 -1 -3];

h1=fliplr(h);%%%%%%%%对序列进行翻转
H=toeplitz(h,[h(1) h1(1:5)]);%%%%%%%%%%%生成循环矩阵
y=H*x';%%%%%%%计算循环卷积序列
%%%%%%%%%%%%FFT
H=fft(h);
X=fft(x);
Y=H.*X;
y1=ifft(Y);

subplot(2,1,1);stem(y);title('直接计算')
subplot(2,1,2);stem(y1);title('DFT计算')



clear all
n1=0:20;
n2=0:10;
h=sinc(0.2*n1);
x=exp(-0.2*n2);
y=conv(x,h);

h1=[h zeros(1,length(x)-1)];
x1=[x zeros(1,length(h)-1)];
H1=fft(h1);
X1=fft(x1);
Y1=H1.*X1;
y1=ifft(Y1);

subplot(2,1,1);stem(y);title('直接计算')
subplot(2,1,2);stem(y1);title('DFT计算')



参考文献

https://www.cnblogs.com/cezorzhao/archive/2013/03/24/2978686.html（用matlab对信号进行傅里叶变换）

 

1 参考文献

[1] 《MATLAB+GUI设计学习手记》;

2 可执行程序生成

花了一天时间浏览了参考文献[1]，初步掌握了matlab gui编程方法。关于gui编程细节不在这里详细介绍，有过matlab编程经验的同学参考《MATLAB+GUI设计学习手记》即可上手。本文档主要介绍matlab环境下standalone应用程序生成方式与没有安装matlab的机器上怎么运行上述程序。步骤如下：

2.1  编译环境

1)  生成standalone应用程序需要mcc编译，首先安装对应编译器。在matlab命令行输入：mbuild –setup； 


2)  输入y； 


3)  选取C语言编译器，推荐“1”为matlab自带的； 


4)  在电脑上添加matlab相关环境变量（我的电脑属性高级环境变量系统变量path），如matlab安装在D盘，则添加“D:\Program Files \MATLAB\R2011b\runtime\win32;D:\Program Files\ MATLAB\R2011b\ bin;”； 

5)  在matlab命令行内输入“mcc –m programfilename.m（后面可以跟若干.m文件，以主函数文件为首，各文件名称以空格分开），即会在当前目录下生成programfilename.exe”；

2.2  无matlab环境下的exe程序运行

在没有安装matlab的机器上运行上述生成的exe程序，按以下步骤即可： 

1)  首先要将\MATLAB701\toolbox\compiler\deploy\win32中的MCRinstaller.exe安装到该计算机上(7.0以前的版本是mglinstaller.exe)； 

2)  其次是要将“MCRinstaller.exe安装目录\runtime\win32”这个路径添加到该计算机的环境变量中，方法同2.1节所述一致；

吐槽一句：运行起来真的很慢！

CUMCM2018 A题解题思路以及代码实现
一些关于CUMCM的思考
CUMCM2018 A题
对于题目的分析
第一问
1.1 总体思路
1.2 模型求解
1.3 Matlab实现
第二问
2.1 总体思路
2.2 模型求解
2.3 Matlab实现
第三问
3.1 总体思路
3.2 模型求解
3.3 Matlab实现
另一种做法
参考文献

一些关于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:N
    M=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:m1
        B(ii)=1+2*r1;A(ii)=-r1;C(ii)=-r1;
        S(ii)=u(ii+1,1);
    end
    for ii=m1+1:m1+m2
        B(ii)=1+2*r2;A(ii)=-r2;C(ii)=-r2;
        S(ii)=u(ii+1,1);
    end
    for ii=m1+m2+1:m1+m2+m3
        B(ii)=1+2*r3;A(ii)=-r3;C(ii)=-r3;
        S(ii)=u(ii+1,1);
    end
    for ii=m1+m2+m3+1:m1+m2+m3+m4-1
        B(ii)=1+2*r4;A(ii)=-r4;C(ii)=-r4;
        S(ii)=u(ii+1,1);
    end
    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);
    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);
    end
    Y=S;
    for k=2:M-1
        Y(k)=S(k)-l(k)*Y(k-1);%中间解
    end
    X=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);
    end
    outputa(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:N
    M=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:m1
        B(ii)=1+2*r1;A(ii)=-r1;C(ii)=-r1;
        S(ii)=u(ii+1,1);
    end
    for ii=m1+1:m1+m2
        B(ii)=1+2*r2;A(ii)=-r2;C(ii)=-r2;
        S(ii)=u(ii+1,1);
    end
    for ii=m1+m2+1:m1+m2+m3
        B(ii)=1+2*r3;A(ii)=-r3;C(ii)=-r3;
        S(ii)=u(ii+1,1);
    end
    for ii=m1+m2+m3+1:m1+m2+m3+m4-1
        B(ii)=1+2*r4;A(ii)=-r4;C(ii)=-r4;
        S(ii)=u(ii+1,1);
    end
        k11=-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);
    end
    Y=S;
    for k=2:M-1
        Y(k)=S(k)-l(k)*Y(k-1);%中间解
    end
    X=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);
    end
    outputa(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:N
    M=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:m1
        B(ii)=1+2*r1;A(ii)=-r1;C(ii)=-r1;
        S(ii)=u(ii+1,1);
    end
    for ii=m1+1:m1+m2
        B(ii)=1+2*r2;A(ii)=-r2;C(ii)=-r2;
        S(ii)=u(ii+1,1);
    end
    for ii=m1+m2+1:m1+m2+m3
        B(ii)=1+2*r3;A(ii)=-r3;C(ii)=-r3;
        S(ii)=u(ii+1,1);
    end
    for ii=m1+m2+m3+1:m1+m2+m3+m4-1
        B(ii)=1+2*r4;A(ii)=-r4;C(ii)=-r4;
        S(ii)=u(ii+1,1);
    end
    k21=-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：追s
        l(k)=A(k)/U(k-1);
        U(k)=B(k)-C(k-1)*l(k);
    end
    Y=S;
    for k=2:M-1
        Y(k)=S(k)-l(k)*Y(k-1);%中间解
    end
    X=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);
    end
    outputa(j,:)=X;
    u(2:M,2)=X; %把结果放入矩阵 u 中
    u(:,1)=u(:,2);% 下一时刻的值
end
if outputa(3300-1,end)>=44;%条件1
    continue;
end
if outputa(:,end)>=47;%条件2
    continue;
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:15
    flag=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:N
    M=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:m1
        B(ii)=1+2*r1;A(ii)=-r1;C(ii)=-r1;
        S(ii)=u(ii+1,1);
    end
    for ii=m1+1:m1+m2
        B(ii)=1+2*r2;A(ii)=-r2;C(ii)=-r2;
        S(ii)=u(ii+1,1);
    end
    for ii=m1+m2+1:m1+m2+m3
        B(ii)=1+2*r3;A(ii)=-r3;C(ii)=-r3;
        S(ii)=u(ii+1,1);
    end
    for ii=m1+m2+m3+1:m1+m2+m3+m4-1
        B(ii)=1+2*r4;A(ii)=-r4;C(ii)=-r4;
        S(ii)=u(ii+1,1);
    end
    k21=-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);
    end
    Y=S;
    for k=2:M-1
        Y(k)=S(k)-l(k)*Y(k-1);%中间解
    end
    X=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);
    end
    outputa(j,:)=X;
    u(2:M,2)=X; %把结果放入矩阵 u 中
    u(:,1)=u(:,2);% 下一时刻值
end
if outputa(1500-1,end)>=44  %条件1
    continue;
end
if outputa(:,end)>=47  %条件2
    continue;
end
num=num+1;
outi(num)=dd;
outj(num)=dd1;
if (outputa(1500-1,end)<44)&(outputa(:,end)<47)
    flag=1;
    break;%得到每个满足要求的第二层厚度对应的第四层最小厚度即退出循环
end
    end
if flag==1
    break;
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