python 一维热传导方程

实战一所学（2018A

• 一维热传导方程
• 1、偏微分方程的解法（pdepe函数）：
• 2、有限差分法

• 写代码遇到的艰难阻碍

• 3、每个时域迭代进行矩阵计算
• matlab较好的取名变量(见名知意）
• 求最大公因数

一维热传导方程

1、偏微分方程的解法（pdepe函数）：

有初值条件和足够的（两个）边值条件
本来想由外层到内层，每层用解偏微分方程的方法解出U(x,t)可惜这题找不到两个边界条件

2、有限差分法

标准的一维热传导方程的解法（有初值条件和足够的（两个）边值条件）

• 类似有限差分法的思路，将时间和空间离散化，

用代码实现以下 推理式

function main_program
clc,clear
format short g
rho=[300,862,74.2,1.18]; %密度
c=[1377,2100,1726,1005]; %比热容
k=[0.082,0.37,0.045,0.028]; %热传导率
x=[0.0006,0.006,0.0036,0.005]; %每层的宽度
% dx=[0.0001,0.001,0.0006,0.001]; %空间步长
gcd = GCD(x); %0.0002
dx=gcd;
dt=0.002; %时间步长
Ten=75;Tbd=37; %初始温度
n=sum(x)/dx+1
m=5400/dt
arr=zeros(m,n);

arr(1,:)=37.*ones(n,1);
arr(:,1)=75.*ones(1,m); %arr(1,1)是75
tmp=dt/(dx*dx);


sj=zeros(size(x));% x是4块的厚度构成的矩阵
sj(1)=x(1);
for i=2:4 
    sj(i)=sj(i-1)+x(i);
end
sj
sj=sj/dx % 3 33 51 76 
%每一块离散化成多少块，再求前缀和，也就是想记录每段接触点的位置
%下面这种方法更好，记录接触点
%上面方法记录每块的长度，还要考虑下标从2开始

% LEN1=int8(x(1)/dx(1))+1; %记录每段接触点的位置
% LEN2=LEN1+x(2)/dx(2);
% LEN3=LEN2+x(3)/dx(3);
% LEN4=LEN3+x(4)/dx(4)

for i=1:m-1
    mu=k(1)/(c(1)*(rho(1)));
    mu=mu*mu;
    for j=2:sj(1)+1
%     arr(i+1,j)= tmp*k(1)*k(1)*(arr(i,j+1)+arr(i,j-1)-2*arr(i,j))/(c(1)*(rho(1)))*(c(1)*(rho(1)))+arr(i,j);
    arr(i+1,j)= tmp*mu*(arr(i,j+1)+arr(i,j-1)-2*arr(i,j))+arr(i,j);
    end
    mu=k(2)/(c(2)*(rho(2)));
    mu=mu*mu;
    for j=sj(1)+2:sj(2)+1  %sj(1)+1  4
        arr(i+1,j)= tmp*mu*(arr(i,j+1)+arr(i,j-1)-2*arr(i,j))+arr(i,j);
    end 
    mu=k(3)/(c(3)*(rho(3)));
    mu=mu*mu;
    for j=sj(2)+2:sj(3)+1
        arr(i+1,j)= tmp*mu*(arr(i,j+1)+arr(i,j-1)-2*arr(i,j))+arr(i,j);
    end 
    mu=k(4)/(c(4)*(rho(4)));
    mu=mu*mu;
    for j=sj(3)+2:sj(4)  %理论上最后一列是sj(4)+1，要按该公式算会数组越界
        arr(i+1,j)= tmp*mu*(arr(i,j+1)+arr(i,j-1)-2*arr(i,j))+arr(i,j);
    end
end

arr(2,:)
arr(3,:)
arr(50,:)
end

很无奈，由于一维热传导方程中$a^2$理解错了，本来就很小的数再加以平方，得出的U(x,t)几乎全部接近37

最后一列全0，因为按照这一个公式算所有的点，最后一列会越界

实际上块与块之间的接触点、热源与x=0处、第四层与皮肤的接触点需要通过别的公式计算，无奈下面三个式子太难推得

写代码遇到的艰难阻碍

1、确定时间步长、空间步长后，将时间、距离热源距离离散为若干个时域和空域（这里时间m份，空间n份）
首先根据边界条件对 离散矩阵的边界进行初始化，t=0，U(x,0)初始化为37，x=0，U(0,t)初始化为75，要注意的是，matlab矩阵下标从1开始，则以空域为例，初始化下标1，1~N对应的矩阵下标是 2 ~ N+1, 故矩阵的大小应初始化为 （m+1，n+1）
2、由于对MATLAB不熟，每个通过计算得到得变量（尤其是矩阵，最好输出看看，比如这里得sj矩阵，一不小心就算出小数，作为数组下标显然会报错位置 2 处的索引无效。数组索引必须为正整数或逻辑值。

3、每个时域迭代进行矩阵计算

超级牛逼的队友推出了某个牛逼式子（我TM直接谢绝观赏），我只负责计算
设如下矩阵为AX=b
对于时域 i ，通过X=A\b，计算出时域 i 的温度，再通过时域 i 的温度 递推出 进行时域 i+1 计算的 b矩阵

首先为得到左侧的A矩阵，主对角线上通过循环，主要是主对角线相邻的两条斜线

本来一开始，初始化上侧，想通过A=A'+A 得到对称矩阵，出错

请注意，通过A=A'+A 得到对称矩阵只适用于 完美对称的矩阵，

像这里，对称的元素所在的行不同，取得$r_i$不同，所以说不是完美对称的，太想当然了

这里就应该分别循环初始化

%文件名称为 skin_T.m
function T = skin_T(Ten,x2,x4)
clc,clear
format short g
rho=[300,862,74.2,1.18]; %密度
c=[1377,2100,1726,1005]; %比热容
k=[0.082,0.37,0.045,0.028]; %热传导率
x=[0.0006,0.006,0.0036,0.005]; %每层的宽度
% dx=[0.0001,0.001,0.0006,0.001]; %空间步长
gcd = GCD(x) %0.0002
dx=0.0002;
dt=0.02; %时间步长
Ten=75;Tbd=37; %初始温度
n=sum(x)/dx % 下标1~ n+1 ,这里不计x=0处的清一色75°
m=5400/dt+1
% U矩阵维度(n,m)

sj=zeros(size(x));
sj(1)=x(1);
for i=2:4 
    sj(i)=sj(i-1)+x(i);
end
sj=sj/dx % 3 33 51 76
% dx
sqrt_a=k./(c.*rho)
r=dt/(dx*dx)*sqrt_a
h=83.6
A=zeros(n,n);
% A=A+-2*eye(n)+diag(ones(1,length(x)-1),1)+diag(ones(1,length(x)-1),-1);

for i=1:n-1
    if(i<=sj(1)+1)A(i,i+1)=-1*r(1);
    elseif(i<=sj(2)+1)A(i,i+1)=-1*r(2);
        elseif(i<=sj(3)+1)A(i,i+1)=-1*r(3);
            elseif(i<=sj(4)+1)A(i,i+1)=-1*r(4);
    end   
end
% A=A'+A ;

for i=2:n
    if(i<=sj(1)+1)A(i,i-1)=-1*r(1);
    elseif(i<=sj(2)+1)A(i,i-1)=-1*r(2);
        elseif(i<=sj(3)+1)A(i,i-1)=-1*r(3);
            elseif(i<=sj(4)+1)A(i,i-1)=-1*r(4);
    end   
end
for i=1:n
    if(i<=sj(1)+1)A(i,i)=1+2*r(1);
    elseif(i<=sj(2)+1)A(i,i)=1+2*r(2);
        elseif(i<=sj(3)+1)A(i,i)=1+2*r(3);
            elseif(i<=sj(4)+1)A(i,i)=1+2*r(4);
    end
end
 A(n,n)=h+k(4)/dx;
 A(n,n-1)=-1*k(4)/dx;
A(1,:);
A(2,:);

 b=zeros(n,1);
 u=zeros(n,m); %m已经算上了t=0的那一列37
%  u=u+37*ones(n,1);
 u(:,1)=37;
 %u(x,t)
% u(1,2)
% u(:,1)


skin=[37];
 for i=2:m %t 270000
     for j=1:n  %x
  if(j==1)b(j)=u(j,i-1)+75*r(1);
  elseif(j==n)
      p=(u(n,i-1)-k(4)/dx*(u(n-1,i-1)-u(n,i-1))/h);
      b(j)=h*p;
      skin=[skin;p];%skin是个列向量
  else b(j)=u(j,i-1);
  end
     end
%      b
 tmp=A\b;

%  u=[u tmp];
%  u=u+tmp ; %u的每一列都会加上tmp列向量
%  u(:,i)=tmp; %正确，但是会自动输出的是整个u矩阵而不是这一列 不利于观察
 u(:,i)=u(:,i)+tmp;
%  for q=1:n
%      u(q,i)=tmp(q);
%  end


u(:,i);
% % i
 end
 disp('--------------------------------');
 u=[u;skin'];

%    xlswrite('work.xlsx',u,'sheet1');
% x=0:dx:sum(x)+dx;
% t=0:dt:5400

 size(u) %n*m
 x=linspace(0,sum(x),n+1) ;
 t=linspace(0,5400,m) ;
 % [X,T]=meshgrid(x,t); 不需要
 mesh(t,x,u); %u(x,t)

shading interp

%  T=skin;
end

%%%%---------------------------子函数如下
%求解最大公因数 Calculate the GCD of an integer array
function result = GCD(array)
n = length(array);
tmp = array;
exp = 0; % Record how many times the array has been multiplied by 10
while ~is_intarray(tmp) % Whether array multiplied by 10's has become integer
    tmp = tmp * 10;
    exp = exp + 1;
end
result = gcd_array(uint32(tmp),n); % Calculate the LCM of the integer array which has been multiplied by 10's
if (exp > 0) % If the array is multiplied by 10's, divide the LCM by 10's in return
    result = double(result) / 10^exp;
end
end
 

function result = gcd_array(array, n)
result = array(1);
for i = 1:n-1
    result = gcd(result,array(i+1)); 
end
end
%判定被乘了10的数组是否已经变成整数了
function result = is_intarray(array)
tmp = round(array);
if (abs(array - tmp) < 10^(-10))
    result = 1;
else
    result = 0;
end
end

要善于输出部分进行测试正误，测试代码如下

%文件名称为 skin_T.m
function T = skin_T(Ten,x2,x4)
% clc,clear
format short g
rho=[300,862,74.2,1.18]; %密度
c=[1377,2100,1726,1005]; %比热容
k=[0.082,0.37,0.045,0.028]; %热传导率
x=[0.0006,0.006,0.0036,0.005]; %每层的宽度
% dx=[0.0001,0.001,0.0006,0.001]; %空间步长
gcd = 0.0002;%GCD(x); %0.0002
dx=0.0002;
dt=0.02; %时间步长
Ten=75;Tbd=37; %初始温度
n=sum(x)/dx ;% 下标1~ n+1 ,这里不计x=0处的清一色75°
m=5400/dt+1;
% U矩阵维度(n,m)
n=5;
sj=zeros(size(x));
sj(1)=x(1);
for i=2:4 
    sj(i)=sj(i-1)+x(i);
end
sj=sj/dx ;% 3 33 51 76
% dx
sqrt_a=k./(c.*rho);
r=dt/(dx*dx)*sqrt_a;
h=83.6;
A=zeros(n,n);


for i=1:n-1
    if(i<=sj(1)+1)A(i,i+1)=-1*r(1);
    elseif(i<=sj(2)+1)A(i,i+1)=-1*r(2);
        elseif(i<=sj(3)+1)A(i,i+1)=-1*r(3);
            elseif(i<=sj(4)+1)A(i,i+1)=-1*r(4);
    end   
end
A=A'+A ;
% 
% for i=2:n
%     if(i<=sj(1)+1)A(i,i-1)=-1*r(1);
%     elseif(i<=sj(2)+1)A(i,i-1)=-1*r(2);
%         elseif(i<=sj(3)+1)A(i,i-1)=-1*r(3);
%             elseif(i<=sj(4)+1)A(i,i-1)=-1*r(4);
%     end   
% end
% for i=1:n
%     if(i<=sj(1)+1)A(i,i)=1+2*r(1);
%     elseif(i<=sj(2)+1)A(i,i)=1+2*r(2);
%         elseif(i<=sj(3)+1)A(i,i)=1+2*r(3);
%             elseif(i<=sj(4)+1)A(i,i)=1+2*r(4);
%     end
% end
%  A(n,n)=h+k(4)/dx;
%  A(n,n-1)=-1*k(4)/dx;
A(1,:);
A(2,:);
A
%  
%  T=skin;
end

matlab较好的取名变量(见名知意）

mu
lamda
theda

求最大公因数

%%%%---------------------------子函数如下
%求解最大公因数 Calculate the GCD of an integer array
function result = GCD(array)
n = length(array);
tmp = array;
exp = 0; % Record how many times the array has been multiplied by 10
while ~is_intarray(tmp) % Whether array multiplied by 10's has become integer
    tmp = tmp * 10;
    exp = exp + 1;
end
result = gcd_array(uint32(tmp),n); % Calculate the LCM of the integer array which has been multiplied by 10's
if (exp > 0) % If the array is multiplied by 10's, divide the LCM by 10's in return
    result = double(result) / 10^exp;
end
end
 

function result = gcd_array(array, n)
result = array(1);
for i = 1:n-1
    result = gcd(result,array(i+1)); 
end
end
%判定被乘了10的数组是否已经变成整数了
function result = is_intarray(array)
tmp = round(array);
if (abs(array - tmp) < 10^(-10))
    result = 1;
else
    result = 0;
end
end