2020国赛A题偏微分方程的数值解法

国赛A题根据Fourier热传导定律和牛顿冷却定律建立了热传导模型，模型的微分方程如下：
$$\{\begin{array}{l}\frac{\partial u}{\partial t}=\frac{k}{c\rho }\frac{{\partial }^{2}u}{\partial x^2}\\ \frac{\partial u}{\partial x}{|}_{x=0}=-\frac{h(T_{炉}-u(0,t))}{k},\\ \frac{\partial u}{\partial x}{|}_{x=l}=-\frac{h(T_{炉}-u(l,t))}{k},\\ u(x,0)=T_0.\end{array}$$
其中， $T_{炉}(x)=$
$$\{\begin{array}{l}{T}_{1-5},\\ \frac{T_6-T_{1-5}}{5}(x-172.5)+T_{1-5},\\ T_6,\\ \frac{T_7-T_6}{5}(x-208)+T_6,\\ T_7,\\ \frac{T_{8-9}-T_7}{5}(x-243.5)+T_7,\\ T_{8-9},\\ \frac{T_{10-11}-T_{8-9}}{5}(x-314.5)+T_{8-9},\\ T_{10-11}\end{array}$$
下面来考虑这个偏微分方程的求解问题，设参数 $\frac{k}{c\rho }$为 $k_1$，参数 $\frac{h}{k}$为 $\frac{1}{k_2}$。

1 向前差分

模型的推导

由一阶泰勒展开和二阶泰勒展开可得：
$$u(x_i,t_{j+1})=u(x_i,t_j)+\frac{\partial u}{\partial t}(x_i,t_j)(t_{j+1}-t_j)+O(t)（1）$$
$$u(x_{i+1},t_j)=u(x_i,t_j)+\frac{\partial u}{\partial x}(x_i,t_j)(x_{i+1}-x_i)+\frac{1}{2}\frac{{\partial }^{2}u}{\partial x^2}(x_i,t_j)(x_{i+1}-x_i{)}^2+O(x^2)（2）$$
$$u(x_{i-1},t_j)=u(x_i,t_j)+\frac{\partial u}{\partial x}(x_i,t_j)(x_{i-1}-x_i)+\frac{1}{2}\frac{{\partial }^{2}u}{\partial x^2}(x_i,t_j)(x_{i-1}-x_i{)}^2+O(x^2)（3）$$
由（1）式可得：
$$\frac{\partial u}{\partial t}(x_i,t_j)=\frac{u(x_i,t_{j+1})-u(x_i,t_j)}{t_{j+1}-t_j}（4）$$
由（2）式和（3）式可得：
$$\frac{{\partial }^{2}u}{\partial x^2}(x_i,t_j)=\frac{u(x_{i+1},t_j)+u(x_{i-1},t_j)-2u(x_i,t_j)}{(x_{i+1}-x_i{)}^2}（5）$$
将(4)式及（5）式代入建立的微分方程中，并将$u(x_i,t_{j+1})$移到式子的一侧，取网格等距划分，可得：
$$u(x_i,t_{j+1})=u(x_i,t_j)+\frac{k}{c\rho }\frac{\Delta t}{(\Delta x{)}^2}[u(x_{i+1},t_j)+u(x_{i-1},t_j)-2u(x_i,t_j)]$$
下面来处理边界条件，这也是用差分法求解偏微分方程关键的一步。此处模型是第三类边界条件，在边界上一阶泰勒展开得：
$$\frac{u(\Delta x,t_j)-u(0,t_j)}{\Delta x}=-\frac{h}{k}(T_{炉}-u(0,t_j))$$
将$u(0,t_j)$移到一边得：
$$u(0,t_j)=u(\Delta x,t_j)+\frac{h\Delta x}{k}(T_{炉}-u(0,t_j))$$
另一边类似处理可得：
$$u(l,t_j)=u(l-\Delta x,t_j)-\frac{h\Delta x}{k}(T_{炉}-u(l,t_j))$$
向前差分法及向后差分法示意图如下：

然后我们就可以编写代码了，写代码的思路就是向前差分的思路：由矩阵中前三个点确定下一个点（由矩阵一列确定下一列），具体代码见下节。

编写代码

代码共分三个板块，其作用在代码注释中已经说的很明白了。

function [TTT,t]=calT2(v,N)
%改编calT1
%计算各个段的炉内温度（分段函数折线）
fdT=[175+273.15 195+273.15 235+273.15 255+273.15 25+273.15];
fdT=fdT-273.15;
% d=0.15;
% x=0:deld:d;%tau--距离向量
%t=0:delt:373; %这个t后面的值不太对,其为题目所给的附件中的时间，可见文件“附件”
t=0:(60*385.5/v)/N:60*385.5/v;
count=zeros(1,9); %各段端点所在位置(温区1-5，温区6，温区7，温区8-9，温区10-11)
for j=1:length(t)
    if t(j)<60*172.5/v && t(j+1)>60*172.5/v
        count(1)=j;
    elseif t(j)<60*177.5/v && t(j+1)>60*177.5/v
        count(2)=j;
    elseif t(j)<60*208/v && t(j+1)>60*208/v
        count(3)=j;
    elseif t(j)<60*213/v && t(j+1)>60*213/v
        count(4)=j;
    elseif t(j)<60*243.5/v && t(j+1)>60*243.5/v
        count(5)=j;    
    elseif t(j)<60*248.5/v && t(j+1)>60*248.5/v
        count(6)=j;
    elseif t(j)<60*314.5/v && t(j+1)>60*314.5/v
        count(7)=j;    
    elseif t(j)<60*319.5/v && t(j+1)>60*319.5/v
        count(8)=j;
    end   
end
%count(9)=length(x); %这里应该是length(t),不过这个值对后续也没什么影响
count(9)=length(t);   
for j=1:length(t)
    if j<=count(1)        
        TTT(j)=fdT(1);   
    elseif j>count(1) && j<count(2)
         %TTT(j)=(fdT(1)-fdT(2))/(t(count(1))-t(count(2)))*t(j)+fdT(1)-t(count(1))*(fdT(1)-fdT(2))/(t(count(1))-t(count(2)));      
         TTT(j)=(fdT(1)+fdT(2))/2;
    elseif j>count(3) && j<count(4)
       % TTT(j)=(fdT(2)-fdT(3))/(t(count(3))-t(count(4)))*t(j)+fdT(2)-t(count(3))*(fdT(2)-fdT(3))/(t(count(3))-t(count(4)));   
       TTT(j)=(fdT(3)+fdT(2))/2;
elseif j>count(5) && j<count(6)
        %TTT(j)=(fdT(3)-fdT(4))/(t(count(5))-t(count(6)))*t(j)+fdT(3)-t(count(5))*(fdT(3)-fdT(4))/(t(count(5))-t(count(6)));  
        TTT(j)=(fdT(3)+fdT(4))/2;
    elseif j>count(7) && j<count(8)
      %  TTT(j)=(fdT(4)-fdT(5))/(t(count(7))-t(count(8)))*t(j)+fdT(4)-t(count(7))*(fdT(4)-fdT(5))/(t(count(7))-t(count(8)));     
      TTT(j)=(fdT(4)+fdT(5))/2;
for j=1:length(t)
    if j<=count(1)        
        TTT(j)=fdT(1);   
    elseif j>count(1) && j<count(2)
         %TTT(j)=(fdT(1)-fdT(2))/(t(count(1))-t(count(2)))*t(j)+fdT(1)-t(count(1))*(fdT(1)-fdT(2))/(t(count(1))-t(count(2)));      
         TTT(j)=(fdT(1)+fdT(2))/2;
    elseif j>count(3) && j<count(4)
       % TTT(j)=(fdT(2)-fdT(3))/(t(count(3))-t(count(4)))*t(j)+fdT(2)-t(count(3))*(fdT(2)-fdT(3))/(t(count(3))-t(count(4)));   
       TTT(j)=(fdT(3)+fdT(2))/2;
    elseif j>count(5) && j<count(6)
        %TTT(j)=(fdT(3)-fdT(4))/(t(count(5))-t(count(6)))*t(j)+fdT(3)-t(count(5))*(fdT(3)-fdT(4))/(t(count(5))-t(count(6)));  
        TTT(j)=(fdT(3)+fdT(4))/2;
    elseif j>count(7) && j<count(8)
      %  TTT(j)=(fdT(4)-fdT(5))/(t(count(7))-t(count(8)))*t(j)+fdT(4)-t(count(7))*(fdT(4)-fdT(5))/(t(count(7))-t(count(8)));     
      TTT(j)=(fdT(4)+fdT(5))/2;
    elseif j>=count(2) && j<=count(3)   
          TTT(j)=fdT(2);        
    elseif j>=count(4) && j<=count(5)           
           TTT(j)=fdT(3);             
    elseif j>=count(6) && j<=count(7)          
           TTT(j)=fdT(4);           
    elseif j>=count(8) 
          TTT(j)=fdT(5);   
    end  
end

function [u,lu,t]=mosol1(v,k1,k2,M,N)
%改编mosol，改成划分格子数的数目来决定划分网格，而不是距离决定
%向前差分计算各个时间点被加热材料上的温度及炉温曲线
%后面优化得到的最佳参数为：
%k1=4.7*10^(-5);k2=1/3900;
%其他参数参考值：deld=0.0005;delt=0.002;v=70;
[T,t]=calT2(v,N);%T--T炉,t--时间向量
% t=0:delt:60*385.5/v;
tau=t(2)-t(1);%时间间隔
d=0.15;
n=200;%n--距离向量分隔密度
x=0:d/M:d;%tau--距离向量
h=x(2)-x(1);%距离间隔
u=zeros(length(x),length(t));
[m,n]=size(u);%m--行数，n--列数
%边界条件
for i=1:m
    u(i,1)=25;
end
for j=1:n-1
    for i=2:m-1
        u(i,j+1)=u(i,j)+(k1*tau)/(h^2)*(u(i-1,j)-2*u(i,j)+u(i+1,j));
    end
    u(1,j+1)=((1/k2)*T(j+1)+u(2,j+1)/h)/(1/h+1/k2);
    u(m,j+1)=(-1/k2*T(j+1)+u(m-1,j+1)/h)/(1/h-1/k2);
end
lu=u(floor(length(x)/2),:);

下面还有一个遍历搜索最优值的代码，语法没错误，就是运行很长时间还是运行不出来，也贴在这里以供参考：

function [kk1,kk2,error]=select()
%优化函数，求得误差最小的参数值
    load shuju tt T;
    v=70;
    kk1=1*10^(-5):0.1*10^(-5):4*10^(-5);
    kk2=1/7000:2*10^(-5):1/3000;
%   kk1=4.7*10^(-5);kk2=1/3900;
    error=zeros(length(kk1),length(kk2));
    delt=0.0005;deld=0.002;
    for i=1:length(kk1)
        for j=1:length(kk2)
            [u1,lu1,t1]=mosol(v,kk1(i),kk2(j),delt,deld);
            nu=interp1(t1,lu1,tt,'spline');
            error(i,j)=norm(nu-T);
        end
    end
    figure
    mesh(kk1,kk2,error);
    xlabel('kk1');
    ylabel('kk2');
    zlabel('error');
end

2 向后差分

模型的推导

向后差分也是由泰勒展开推导得的，与向前差分的区别在于它需要求解一个线性方程组才能求能由矩阵中前面已知的一列推得后一列的值。
$$\frac{u(x_i,t_j)-u(x_i,t_{j-1})}{\Delta t}=\frac{k}{c\rho (\Delta x{)}^2}[u(x_{i+1},t_j)+u(x_{i-1},t_j)-2u(x_i,t_j)]$$整理得：
$$u(x_i,t_{j-1})=(1+2r)u(x_i,t_j)-r(u(x_{i+1},t_j)+u(x_{i-1},u(x_{i-1},t_j)))\text{其中，}r=\frac{k}{c\rho }\frac{\Delta t}{(\Delta x{)}^2}$$
据此，我们可建立线性方程组，但仍需处理边界条件。
$$\left(\begin{array}{c}1+2r-r0\cdots 0\\ -r1+2r-r\ddots ⋮\\ 0-r1+2r\ddots 0\\ ⋮\ddots \ddots \ddots -r\\ 0\cdots 0-r1+2r\end{array}\right)\left(\begin{array}{c}u(1,j)\\ ⋮\\ u(M-1,j)\end{array}\right)=\left(\begin{array}{c}u(1,j-1)\\ ⋮\\ u(M-1,j-1)\end{array}\right)+r\left(\begin{array}{c}u(0,j)\\ 0\\ ⋮\\ 0\\ u(M,j)\end{array}\right)$$
其中，$u(0,j)$和$u(M,j)$在第三类边界条件下是未知的（在第一类边界条件下是已知的），但可根据边界条件确定其与领接点的关系式，由上面向前差分法的边界条件可得：
$$u(0,j)=\frac{u(1,j)+\frac{h}{k}\Delta x{T}_{炉}}{1+\frac{h}{k}\Delta x}$$
$$u(M,j)=\frac{u(M,j)-\frac{h}{k}\Delta x{T}_{炉}}{1-\frac{h}{k}\Delta x}$$
将上式中$u(1,j)$,$u(M,j)$移到线性方程组左边整理得：
$$\left(\begin{array}{c}1+r-r0\cdots 0\\ -r1+2r-r\ddots ⋮\\ 0-r1+2r\ddots 0\\ ⋮\ddots \ddots \ddots -r\\ 0\cdots 0-r1+r\end{array}\right)\left(\begin{array}{c}u(1,j)\\ ⋮\\ u(M-1,j)\end{array}\right)=\left(\begin{array}{c}u(1,j-1)\\ ⋮\\ u(M-1,j-1)\end{array}\right)+r\left(\begin{array}{c}\frac{\frac{h}{k}\Delta x{T}_{炉}}{1+\frac{h}{k}\Delta x}\\ 0\\ ⋮\\ 0\\ -\frac{\frac{h}{k}\Delta x{T}_{炉}}{1-\frac{h}{k}\Delta x}\end{array}\right)$$
这样就得到前面一列的温度值和后面一列温度值的关系，这是一个三对角占优方程组，可用追赶法（LU分解）求解，要依此求解，只需放在循环环境下即可。

编写代码

具体代码如下：

  function [u,lu,t]=mosol_xhcf1(k1,k2,v,M,N)
 %改编xhxf
 %向后差分法求数值解,与之相关的代码为calT1(计算各个段的炉内温度（分段函数折线）),zhong_tridiag（追赶法求解三对角线性系统,a,b,c对应三对角矩阵上的元素，f是线性系统的右边项）
 %参数：k1=4.7*10^(-5);k2=1/3900;N=6600;M=300;v=70;
 %边界条件
 %u(z,0)=25 初值条件
 [TTT,t]=calT2(v,N);
 d=0.15; %单位为cm
 delt=(60*385.5/v)/N;deld=d/M; %M=100;N=5000;%M,N涉及到了差分格式的敏感性，根据具体的问题来进行相应的调整
 x=0:deld:d;%t=0:delt:385.5/v;
 u=zeros(length(x),length(t));
 %M=length(x)-1;N=length(t)-1;
 %初值
 for i=1:M+1
 u(i,1)=25;
 end
 %% 已经有追赶法的程序，只需要把输入整理进去
 r=k1*delt/(deld^2);
 aa(1:M-2)=-r;
 cc=aa;
 bb(1)=1+r;bb(2:M-2)=1+2*r;bb(M-1)=1+r;
 for j=2:N+1
 dd=[u(2,j-1)+r*(1/k2)*deld*TTT(j)/(1+(1/k2)*deld);u(3:M-1,j-1);
 u(M,j-1)-r*deld*(1/k2)*TTT(j)/(1+(1/k2)*deld)];
 u(2:M,j)=zhong_tridiag(aa,bb,cc,dd);
 end
 lu=u(floor(M/2),:);
 figure
 plot(t,lu)

function x =zhong_tridiag(a,b,c,f)
%追赶法求解三对角线性系统,a,b,c对应三对角矩阵上的元素，f是线性系统的右边项
n=length(f);
beta=zeros(n-1,1); y=zeros(n,1); x=zeros(n,1);%初始定义
beta(1)=c(1)/b(1);
for j=2:n-1
beta(j)=c(j)/(b(j)-a(j-1)*beta(j-1));%对应step 1
end
y(1)=f(1)/b(1);
for j=2:n
y(j)=(f(j)-a(j-1)*y(j-1))/(b(j)-a(j-1)*beta(j-1));%对应step 2
end
x(n)=y(n);
for j=n-1:-1:1
x(j)=y(j)-beta(j)*x(j+1);%对应step 3
end 
end

其中用到的calT2和前面的向前差分法的代码一样。

3 两种方法比较

可以看到向后差分法画出的效果图基本不变，而向前差分法画出的效果图在N较小的时候效果很差，这也说明了向后差分法比向前差分法更稳定。 此处求解得到的解与题目中所给的数据之间还存在较大的差距（此处向后差分差距更大，可能是向后差分没有调参数的问题），可能是没有调好参数的问题（可以试试每个温度段都设置一个参数，这样得到的效果应该会好很多），但与实验数据差不多从一个方面说明了我们求解偏微分方程的数值解法的代码是正确的。
画图及运行前面编写的代码的代码如下：

function oper()
k1=4.7*10^(-5);k2=1/3900;N=66000;M=300;v=70;
[u,lu,t]=mosol_xhcf1(k1,k2,v,M,N);
plot(t,lu)
hold on
[u,lu,t]=mosol1(v,k1,k2,M,N);
figure
plot(t,lu)
N=6600;M=300;
[u,lu,t]=mosol_xhcf1(k1,k2,v,M,N);
figure(1)
plot(t,lu)
figure(2)
[u,lu,t]=mosol1(v,k1,k2,M,N);
plot(t,lu)
legend('M=300,N=66000向前差分','M=300,N=6600向前差分')
N=166000;M=300;
[u,lu,t]=mosol_xhcf1(k1,k2,v,M,N);
figure(1)
plot(t,lu)
legend('M=300,N=66000向后差分','M=300,N=6600向后差分','M=300,N=166000向后差分')
figure
[u,lu,t]=mosol1(v,k1,k2,M,N);
plot(t,lu)
legend('M=300,N=166000向前差分')

附录

草稿纸(逻辑推理部分)