-
題目
在高溫環境下工作時,人們需要穿着專用服裝以避免灼傷。專用服裝通常由三層織物材料構成,記為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. 假人皮膚外側的測量溫度
題目下載下傳位址:http://special.univs.cn/service/jianmo/sxjmtmhb/2018/0709/1179372.shtml
- 解答如下:
-
假設:
假設 1: 假設在整個過程中外部溫度保持不變;
假設 2: 假設在整個的過程中隻考慮熱傳導這一種熱傳遞方式, 熱對流以及熱輻射可以
忽略;
假設 3: 假設織物材料無褶皺且其表面溫度處處相同;
假設 4: 假設在實驗中熱傳導的方向與假人人體垂直;
假設 5: 假設假人可以承受實驗中所達到的所有可能溫度;
假設 6: 假設每一層的參數固定不變, 即不會随着溫度的改化發生變化;
假設 7: 假設專用服裝的初始溫度與假人體内溫度一樣, 為 37℃;
問題(1)
為了求出溫度的分布, 我們需要通過建立溫度——外部到假人皮膚之間的距離——時間的數學模型(示意圖如下圖)
用微元法來建立熱防護服系統的熱量傳遞模型
微元[x,x+Δx]在時間段[t,t+Δt]内吸收的熱量為:
Fourier熱傳導定律 Q1=Q2上述方程在各層之間的區間内成立
邊界條件:
左邊界
右邊界 交界面 熱防護服系統熱量傳遞的偏微分方程模型:再通過有限差分格式求解:
ks為第四層與皮膚之間的熱交換系數 ks=8.36
ke為第I層與高溫環境之間的熱交換系數 ke=117.40
(2)這裡用的二分法,在滿足條件的情況下,求最小的厚度
L2=17.6609mm
(3)這問要求第二層,第四層的最優厚度,雙目标最優函數,看了一下各路大神的思路,無非就是遺傳算法,暴力窮舉,個人認為這裡将主要考慮的是産品的成本,皮膚外側溫度随着第二層、第四層厚度的增加而降低的,這個是可以驗證的,第四層越厚,第二層必然就會薄一些, 成本自然就低一些,而第四層厚度的增加,所需要的服裝材料幾乎沒有變化。是以我這裡将第四層的厚度取最大,将第四層厚度确定之後,這個問題就轉化成第二問的問題了,這樣就很簡單了。
L2=19.2813mm
L4=6.4mm
- 具體代碼如下:
問題(1)
%問題一的主函數 T1.m
clc
clear;
T=1; %時間步長
m=[6 60 36 50];%将每一層細分為0.1mm每塊
%m=[30 300 180 250] %将每一層細分為0.02mm每塊
ms=[0 0 0 0]; % ms表示總的分層數
ms(1,1)=m(1,1);
for j=2:4
ms(1,j)=ms(1,j-1)+m(1,j);
end
l=[0.6 6.6 10.2 15.2]*10^(-3); %l1 l2 l3 l4z的坐标
p=[300 862 74.2 1.18]; %各層密度
c=[1377 2100 1726 1005]; %比熱
k=[0.082 0.37 0.045 0.028]; %熱傳導率
u0=37.0; %初始溫度
ue=75;
%求空間步長
x=zeros(1,4);
x(1,1)=l(1,1)/m(1,1);
for j=2:4
x(1,j)=(l(1,j)-l(1,j-1))/m(1,j);
end
% M(1,j) 從第1層到第j層的等分數
M(1,1)=m(1,1);
for j=2:4
M(1,j)=M(1,j-1)+m(1,j);
end
%求式(10)裡的λ,用b表示λ
for j=1:4
b(1,j)=(k(1,j)*T)/(c(1,j)*p(1,j)*(x(1,j)^2));
end
ke=117.4035;
e=(ke*x(1,1)/k(1,1)); %求式(12)的ue 用e表示
ks=8.36;
us=(ks*x(1,4)/k(1,4)); %求式(13)的us
for j=1:4
v(1,j)=k(1,j)/x(1,j); %求式(14)的v
end
u=zeros(5400/T+1,ms(1,4)+1); %建立一個以距離為x軸,時間為y軸的坐标軸
u(1,:)=37; %時間t=0時,各點溫度為37度
f=zeros(ms(1,4)+1,1); %AX=f
x01=e+1; %l0的邊界條件
x02=-1;
x1=[-b(1,1) 1+2*b(1,1) -b(1,1)]; %第一層
x2=[-b(1,2) 1+2*b(1,2) -b(1,2)]; %第二層
x3=[-b(1,3) 1+2*b(1,3) -b(1,3)]; %第三層
x4=[-b(1,4) 1+2*b(1,4) -b(1,4)]; %第四層
x12=[-v(1,1) v(1,1)+v(1,2) -v(1,2)]; %12層邊界
x23=[-v(1,2) v(1,2)+v(1,3) -v(1,3)];%23層邊界
x34=[-v(1,3) v(1,3)+v(1,4) -v(1,4)];%34層邊界
x51=-1; %l4的邊界條件
x52=1+us;
A=zeros(ms(1,4)+1,ms(1,4)+1);
A(1,1)=x01;
A(1,2)=x02;
f(1,1)=e*ue;
for i=2:ms(1,1)
A(i,i-1:i+1)=x1;
f(i,1)=37;
end
A(ms(1,1)+1,ms(1,1):ms(1,1)+2)=x12;
f(ms(1,1)+1,1)=0;
for i=ms(1,1)+2:ms(1,2)
A(i,i-1:i+1)=x2;
f(i,1)=37;
end
A(ms(1,2)+1,ms(1,2):ms(1,2)+2)=x23;
f(ms(1,2)+1,1)=0;
for i=ms(1,2)+2:ms(1,3)
A(i,i-1:i+1)=x3;
f(i,1)=37;
end
A(ms(1,3)+1,ms(1,3):ms(1,3)+2)=x34;
f(ms(1,3)+1,1)=0;
for i=ms(1,3)+2:ms(1,4)
A(i,i-1:i+1)=x4;
f(i,1)=37;
end
A(ms(1,4)+1,ms(1,4))=x51;
A(ms(1,4)+1,ms(1,4)+1)=x52;
f(ms(1,4)+1,1)=us*u0;
%追趕法
a=zeros(1,ms(1,4)+1);%下對角元素
a(1,1)=0;
for i=2:ms(1,4)+1
a(1,i)=A(i,i-1);
end
c=zeros(1,ms(1,4)+1);%上對角元素
c(1,ms(1,4)+1)=0;
for i=1:ms(1,4)
c(1,i)=A(i,i+1);
end
b=zeros(1,ms(1,4)+1);
for i=1:ms(1,4)+1
b(1,i)=A(i,i);
end
y=chase(a,b,c,f);
ans=zeros(5400/T+1,ms(1,4)+1); %結果矩陣
ans(1,:)=u0;
ans(2,:)=y;
%循環求結果
for i=2:5400/T
f=y';
f(1,1)=e*ue;
f(ms(1,1)+1,1)=0;
f(ms(1,2)+1,1)=0;
f(ms(1,3)+1,1)=0;
f(ms(1,4)+1,1)=us*u0;
y=chase(a,b,c,f);
ans(i+1,:)=y;
end
xlswrite('problem1.xlsx',ans,'Sheet2','A1');
figure(1)
%溫度随時間、空間的分布圖
t=0:T:5400;
y=0:ms(1,4);
[b1,a1]=meshgrid(y,t);
mesh(b1,a1,ans);
title('溫度随時間、空間的分布圖');
%xlabel('空間/0.02mm');
xlabel('空間/0.1mm');
ylabel('時間/s');
zlabel('溫度/℃');
figure(2)
%分界面的溫度随時間的變化圖
t=0:T:5400;
%左邊界
l1=ans(:,1);
%第一層和第二層邊界
l2=ans(:,ms(1,1)+1);
%第二層和第三層邊界
l3=ans(:,ms(1,2)+1);
%第三層和第四層邊界
l4=ans(:,ms(1,3)+1);
%右邊界I、II、III、IV
l5=ans(:,ms(1,4)+1);
plot(t,l1,'m-',t,l2,'g-',t,l3,'r-',t,l4,'c-',t,l5,'k-');
xlabel('時間/s');
ylabel('溫度/℃');
title('分界面的溫度随時間的變化');
legend('左邊界','I、II層邊界','II、III層邊界','III、IV層邊界','右邊界');
%追趕法 chase.m
function x=chase(a,b,c,f)
%求解線性方程組 Ax=f, 其中 A是三對角陣
%a是矩陣 A的下對角線元素 a(1)=0
%b是矩陣 A的對角線元素
%c是矩陣 A的上對角線元素 c(n)=0
%f是方程組的右端向量
n=length(f);
x=zeros(1,n);
y=zeros(1,n);
d=zeros(1,n);
z= zeros(1,n);
%預處理
d(1)=b(1);
for i=1:n-1
z(i)=c(i)/d(i);
d(i+1)=b(i+1)-a(i+1)*z(i);
end
%追的過程
y(1)=f(1)/d(1);
for i=2:n
y(i)=(f(i)-a(i)*y(i-1))/d(i);
end
%趕的過程
x(n)=y(n);
for i=n-1:-1:1
x(i)=y(i)-z(i)*x(i+1);
end
% 求ke,ks quedingxishu.m
clc,clear;
a=load('TT.txt'); %TT.txt 存儲的是附件2的所有數字資訊
ue=75; %環境溫度
us=37; %皮膚内測溫度
ks=8.36; %與皮膚的熱交換系數
keq=1;%最合适的 與環境的熱交換系數
ksq=1;%最合适的 與皮膚的熱交換系數
min=100000; %初始值
k=[0.082 0.37 0.045 0.028]; %熱傳導率
l=[0.6 6.6 10.2 15.2]*10^-3;
sum=l(1,1)/k(1,1);
for j=2:4
sum=sum+(l(1,j)-l(1,j-1))/k(1,j);
end
um=48.08; %穩态溫度
time=a(:,2);
time=time';
for ks=7.8:0.02:8.8
ke=1/((ue-um)/(ks*(um-us))-sum); %與外界的熱交換系數
w=ul4(ke,ks);
min1=0;
for i=1:5041
min1=min1+1/2*(time(1,i)-w(1,i))^2;
end
if min1<min
min=min1;
keq=ke;
ksq=ks;
end
end
min
keq %最合适的 與環境的熱交換系數
ksq %最合适的 與皮膚的熱交換系數
%問題1 傳回穩态的溫度用于求ke,ks的值 ul4.m
%問題1 傳回穩态的溫度用于求ke,ks ul4.m
function [q] = ul4( ke,ks )
T=1; %時間步長
m=[6 60 36 50];%将每一層細分為0.1mm每塊
%m=[30 300 180 250] %将每一層細分為0.02mm每塊
ms=[0 0 0 0]; % ms表示總的分層數
ms(1,1)=m(1,1);
for j=2:4
ms(1,j)=ms(1,j-1)+m(1,j);
end
l=[0.6 6.6 10.2 15.2]*10^(-3); %l1 l2 l3 l4z的坐标
p=[300 862 74.2 1.18]; %各層密度
c=[1377 2100 1726 1005]; %比熱
k=[0.082 0.37 0.045 0.028]; %熱傳導率
u0=37.0; %初始溫度
ue=75;
%求空間步長
x=zeros(1,4);
x(1,1)=l(1,1)/m(1,1);
for j=2:4
x(1,j)=(l(1,j)-l(1,j-1))/m(1,j);
end
% M(1,j) 從第1層到第j層的等分數
M(1,1)=m(1,1);
for j=2:4
M(1,j)=M(1,j-1)+m(1,j);
end
%求式(10)裡的λ,用b表示λ
for j=1:4
b(1,j)=(k(1,j)*T)/(c(1,j)*p(1,j)*(x(1,j)^2));
end
e=(ke*x(1,1)/k(1,1)); %求式(12)的ue 用e表示
us=(ks*x(1,4)/k(1,4)); %求式(13)的us
for j=1:4
v(1,j)=k(1,j)/x(1,j); %求式(14)的v
end
u=zeros(5400/T+1,ms(1,4)+1); %建立一個以距離為x軸,時間為y軸的坐标軸
u(1,:)=37; %時間t=0時,各點溫度為37度
f=zeros(ms(1,4)+1,1); %AX=f
x01=e+1; %l0的邊界條件
x02=-1;
x1=[-b(1,1) 1+2*b(1,1) -b(1,1)]; %第一層
x2=[-b(1,2) 1+2*b(1,2) -b(1,2)]; %第二層
x3=[-b(1,3) 1+2*b(1,3) -b(1,3)]; %第三層
x4=[-b(1,4) 1+2*b(1,4) -b(1,4)]; %第四層
x12=[-v(1,1) v(1,1)+v(1,2) -v(1,2)]; %12層邊界
x23=[-v(1,2) v(1,2)+v(1,3) -v(1,3)];%23層邊界
x34=[-v(1,3) v(1,3)+v(1,4) -v(1,4)];%34層邊界
x51=-1; %l4的邊界條件
x52=1+us;
A=zeros(ms(1,4)+1,ms(1,4)+1);
A(1,1)=x01;
A(1,2)=x02;
f(1,1)=e*ue;
for i=2:ms(1,1)
A(i,i-1:i+1)=x1;
f(i,1)=37;
end
A(ms(1,1)+1,ms(1,1):ms(1,1)+2)=x12;
f(ms(1,1)+1,1)=0;
for i=ms(1,1)+2:ms(1,2)
A(i,i-1:i+1)=x2;
f(i,1)=37;
end
A(ms(1,2)+1,ms(1,2):ms(1,2)+2)=x23;
f(ms(1,2)+1,1)=0;
for i=ms(1,2)+2:ms(1,3)
A(i,i-1:i+1)=x3;
f(i,1)=37;
end
A(ms(1,3)+1,ms(1,3):ms(1,3)+2)=x34;
f(ms(1,3)+1,1)=0;
for i=ms(1,3)+2:ms(1,4)
A(i,i-1:i+1)=x4;
f(i,1)=37;
end
A(ms(1,4)+1,ms(1,4))=x51;
A(ms(1,4)+1,ms(1,4)+1)=x52;
f(ms(1,4)+1,1)=us*u0;
%追趕法
a=zeros(1,ms(1,4)+1);%下對角元素
a(1,1)=0;
for i=2:ms(1,4)+1
a(1,i)=A(i,i-1);
end
c=zeros(1,ms(1,4)+1);%上對角元素
c(1,ms(1,4)+1)=0;
for i=1:ms(1,4)
c(1,i)=A(i,i+1);
end
b=zeros(1,ms(1,4)+1);
for i=1:ms(1,4)+1
b(1,i)=A(i,i);
end
y=chase(a,b,c,f);
ans=zeros(5400/T+1,ms(1,4)+1); %結果矩陣
ans(1,:)=u0;
ans(2,:)=y;
%循環求結果
for i=2:5400/T
f=y';
f(1,1)=e*ue;
f(ms(1,1)+1,1)=0;
f(ms(1,2)+1,1)=0;
f(ms(1,3)+1,1)=0;
f(ms(1,4)+1,1)=us*u0;
y=chase(a,b,c,f);
ans(i+1,:)=y;
end
q=ans(:,ms(1,4)+1);
q=q';
end
問題(2)
%問題2的主函數 T2.m
clc,clear;
len=0.1; %二分法結束的條件
min=0.6;
max=25;
ans=0;%結果
q=fun2(min);%求解最小厚度情況下是否滿足限制條件
if q(1,3601)<=47&&q(1,3301)<=44
ans=min;
fprintf('最小厚度情況下滿足限制條件');
end
q=fun2(max);%求解最大厚度情況下是否滿足限制條件
if q(1,3601)>47||q(1,3301)>44
fprintf('最大厚度情況下不滿足限制條件,此題無解');
end
while(max-min>len)
mid=(min+max)/2;
q=fun2(mid);
if q(1,3601)<=47&&q(1,3301)<=44
max=mid;
else
min=mid;
end
end
ans=max
%對于不同的L2,傳回皮膚外側溫度變化 fun2.m
function [q]=fun2(l2)
%l2 傳入II層的厚度
T=1;%時間步長
s=ceil(l2/0.1);
m=zeros(1,4);
%将每一層細分為0.1mm每塊
m(1,1)=6;
m(1,2)=s;
m(1,3)=36;
m(1,4)=55;
ms=[0 0 0 0]; % ms表示總的分層數
ms(1,1)=m(1,1);
for j=2:4
ms(1,j)=ms(1,j-1)+m(1,j);
end
l=[0.6 l2+0.6 l2+4.2 l2+9.7]*10^(-3); %l1 l2 l3 l4的坐标
p=[300 862 74.2 1.18]; %各層密度
c=[1377 2100 1726 1005]; %比熱
k=[0.082 0.37 0.045 0.028]; %熱傳導率
u0=37.0; %初始溫度
ue=65; %環境溫度
%求空間步長
x=zeros(1,4);
x(1,1)=l(1,1)/m(1,1);
for j=2:4
x(1,j)=(l(1,j)-l(1,j-1))/m(1,j);
end
% M(1,j) 從第1層到第j層的等分數
M(1,1)=m(1,1);
for j=2:4
M(1,j)=M(1,j-1)+m(1,j);
end
%求式(10)裡的λ,用b表示λ
for j=1:4
b(1,j)=(k(1,j)*T)/(c(1,j)*p(1,j)*(x(1,j)^2));
end
ke=117.4035;
e=(ke*x(1,1)/k(1,1)); %求式(12)的ue 用e表示
ks=8.36;
us=(ks*x(1,4)/k(1,4)); %求式(13)的us
for j=1:4
v(1,j)=k(1,j)/x(1,j); %求式(14)的v
end
u=zeros(5400/T+1,ms(1,4)+1); %建立一個以距離為x軸,時間為y軸的坐标軸
u(1,:)=37; %時間t=0時,各點溫度為37度
f=zeros(ms(1,4)+1,1); %AX=f
x01=e+1; %l0的邊界條件
x02=-1;
x1=[-b(1,1) 1+2*b(1,1) -b(1,1)]; %第一層
x2=[-b(1,2) 1+2*b(1,2) -b(1,2)]; %第二層
x3=[-b(1,3) 1+2*b(1,3) -b(1,3)]; %第三層
x4=[-b(1,4) 1+2*b(1,4) -b(1,4)]; %第四層
x12=[-v(1,1) v(1,1)+v(1,2) -v(1,2)]; %12層邊界
x23=[-v(1,2) v(1,2)+v(1,3) -v(1,3)];%23層邊界
x34=[-v(1,3) v(1,3)+v(1,4) -v(1,4)];%34層邊界
x51=-1; %l4的邊界條件
x52=1+us;
A=zeros(ms(1,4)+1,ms(1,4)+1);
A(1,1)=x01;
A(1,2)=x02;
f(1,1)=e*ue;
for i=2:ms(1,1)
A(i,i-1:i+1)=x1;
f(i,1)=37;
end
A(ms(1,1)+1,ms(1,1):ms(1,1)+2)=x12;
f(ms(1,1)+1,1)=0;
for i=ms(1,1)+2:ms(1,2)
A(i,i-1:i+1)=x2;
f(i,1)=37;
end
A(ms(1,2)+1,ms(1,2):ms(1,2)+2)=x23;
f(ms(1,2)+1,1)=0;
for i=ms(1,2)+2:ms(1,3)
A(i,i-1:i+1)=x3;
f(i,1)=37;
end
A(ms(1,3)+1,ms(1,3):ms(1,3)+2)=x34;
f(ms(1,3)+1,1)=0;
for i=ms(1,3)+2:ms(1,4)
A(i,i-1:i+1)=x4;
f(i,1)=37;
end
A(ms(1,4)+1,ms(1,4))=x51;
A(ms(1,4)+1,ms(1,4)+1)=x52;
f(ms(1,4)+1,1)=us*u0;
%追趕法
a=zeros(1,ms(1,4)+1);%下對角元素
a(1,1)=0;
for i=2:ms(1,4)+1
a(1,i)=A(i,i-1);
end
c=zeros(1,ms(1,4)+1);%上對角元素
c(1,ms(1,4)+1)=0;
for i=1:ms(1,4)
c(1,i)=A(i,i+1);
end
b=zeros(1,ms(1,4)+1);
for i=1:ms(1,4)+1
b(1,i)=A(i,i);
end
y=chase(a,b,c,f);
ans=zeros(5400/T+1,ms(1,4)+1); %結果矩陣
ans(1,:)=u0;
ans(2,:)=y;
%循環求結果
for i=2:5400/T
f=y';
f(1,1)=e*ue;
f(ms(1,1)+1,1)=0;
f(ms(1,2)+1,1)=0;
f(ms(1,3)+1,1)=0;
f(ms(1,4)+1,1)=us*u0;
y=chase(a,b,c,f);
ans(i+1,:)=y;
end
q=ans(:,ms(1,4)+1);
q=q';
end
%求假人皮膚外側穩定溫度随II層厚度的變化 draw.m
clc,clear;
i=0;
for l2=0.6:0.1:25
q=fun2(l2);
i=i+1;
a(1,i)=q(1,length(q));
end
l2=0.6:0.1:25;
plot(l2,a,'c');
title('假人皮膚外側穩定溫度随II層厚度的變化');
xlabel('II層厚度/mm');
ylabel('假人皮膚外側穩定溫度/攝氏度');
問題(3)
%問題3的主函數 T3.m
clc,clear;
len=0.1; %二分法結束的條件
min=0.6;
max=25;
ans=0;%結果
q=fun3(min);%求解最小厚度情況下是否滿足限制條件
if q(1,1801)<=47&&q(1,1501)<=44
ans=min;
fprintf('最小厚度情況下滿足限制條件');
end
q=fun3(max);%求解最大厚度情況下是否滿足限制條件
if q(1,1801)>47||q(1,1501)>44
fprintf('最大厚度情況下不滿足限制條件,此題無解');
end
while(max-min>len)
mid=(min+max)/2;
q=fun3(mid);
if q(1,1801)<=47&&q(1,1501)<=44
max=mid;
else
min=mid;
end
end
ans=max
%傳入不同的L2,傳回皮膚外側的溫度變化矩陣 fun3.m
function [q]=fun3(l2)
%l2 傳入II層的厚度
T=1;%時間步長
s=ceil(l2/0.1);
m=zeros(1,4);
%将每一層細分為0.1mm每塊
m(1,1)=6;
m(1,2)=s;
m(1,3)=36;
m(1,4)=64; %第四層取6.4mm
ms=[0 0 0 0]; % ms表示總的分層數
ms(1,1)=m(1,1);
for j=2:4
ms(1,j)=ms(1,j-1)+m(1,j);
end
l=[0.6 l2+0.6 l2+4.2 l2+10.6]*10^(-3); %l1 l2 l3 l4的坐标
p=[300 862 74.2 1.18]; %各層密度
c=[1377 2100 1726 1005]; %比熱
k=[0.082 0.37 0.045 0.028]; %熱傳導率
u0=37.0; %初始溫度
ue=80; %環境溫度
%求空間步長
x=zeros(1,4);
x(1,1)=l(1,1)/m(1,1);
for j=2:4
x(1,j)=(l(1,j)-l(1,j-1))/m(1,j);
end
% M(1,j) 從第1層到第j層的等分數
M(1,1)=m(1,1);
for j=2:4
M(1,j)=M(1,j-1)+m(1,j);
end
%求式(10)裡的λ,用b表示λ
for j=1:4
b(1,j)=(k(1,j)*T)/(c(1,j)*p(1,j)*(x(1,j)^2));
end
ke=117.4035;
e=(ke*x(1,1)/k(1,1)); %求式(12)的ue 用e表示
ks=8.36;
us=(ks*x(1,4)/k(1,4)); %求式(13)的us
for j=1:4
v(1,j)=k(1,j)/x(1,j); %求式(14)的v
end
u=zeros(5400/T+1,ms(1,4)+1); %建立一個以距離為x軸,時間為y軸的坐标軸
u(1,:)=37; %時間t=0時,各點溫度為37度
f=zeros(ms(1,4)+1,1); %AX=f
x01=e+1; %l0的邊界條件
x02=-1;
x1=[-b(1,1) 1+2*b(1,1) -b(1,1)]; %第一層
x2=[-b(1,2) 1+2*b(1,2) -b(1,2)]; %第二層
x3=[-b(1,3) 1+2*b(1,3) -b(1,3)]; %第三層
x4=[-b(1,4) 1+2*b(1,4) -b(1,4)]; %第四層
x12=[-v(1,1) v(1,1)+v(1,2) -v(1,2)]; %12層邊界
x23=[-v(1,2) v(1,2)+v(1,3) -v(1,3)];%23層邊界
x34=[-v(1,3) v(1,3)+v(1,4) -v(1,4)];%34層邊界
x51=-1; %l4的邊界條件
x52=1+us;
A=zeros(ms(1,4)+1,ms(1,4)+1);
A(1,1)=x01;
A(1,2)=x02;
f(1,1)=e*ue;
for i=2:ms(1,1)
A(i,i-1:i+1)=x1;
f(i,1)=37;
end
A(ms(1,1)+1,ms(1,1):ms(1,1)+2)=x12;
f(ms(1,1)+1,1)=0;
for i=ms(1,1)+2:ms(1,2)
A(i,i-1:i+1)=x2;
f(i,1)=37;
end
A(ms(1,2)+1,ms(1,2):ms(1,2)+2)=x23;
f(ms(1,2)+1,1)=0;
for i=ms(1,2)+2:ms(1,3)
A(i,i-1:i+1)=x3;
f(i,1)=37;
end
A(ms(1,3)+1,ms(1,3):ms(1,3)+2)=x34;
f(ms(1,3)+1,1)=0;
for i=ms(1,3)+2:ms(1,4)
A(i,i-1:i+1)=x4;
f(i,1)=37;
end
A(ms(1,4)+1,ms(1,4))=x51;
A(ms(1,4)+1,ms(1,4)+1)=x52;
f(ms(1,4)+1,1)=us*u0;
%追趕法
a=zeros(1,ms(1,4)+1);%下對角元素
a(1,1)=0;
for i=2:ms(1,4)+1
a(1,i)=A(i,i-1);
end
c=zeros(1,ms(1,4)+1);%上對角元素
c(1,ms(1,4)+1)=0;
for i=1:ms(1,4)
c(1,i)=A(i,i+1);
end
b=zeros(1,ms(1,4)+1);
for i=1:ms(1,4)+1
b(1,i)=A(i,i);
end
y=chase(a,b,c,f);
ans=zeros(5400/T+1,ms(1,4)+1); %結果矩陣
ans(1,:)=u0;
ans(2,:)=y;
%循環求結果
for i=2:5400/T
f=y';
f(1,1)=e*ue;
f(ms(1,1)+1,1)=0;
f(ms(1,2)+1,1)=0;
f(ms(1,3)+1,1)=0;
f(ms(1,4)+1,1)=us*u0;
y=chase(a,b,c,f);
ans(i+1,:)=y;
end
q=ans(:,ms(1,4)+1);
q=q';
end
代碼都是自己花時間碼出來的,問題應該是沒有的,不足就是重複性比較高。非常感謝大家的關注與提問,因本題是2018國賽的題目,網上的優秀論文也比較多,花個幾天時間去了解,算法也比較容易實作。因為作者在準備考研,時間比較緊張,不能及時回複大家,還請多多見諒。
作者有幸參加了2019年高教社杯數學模組化,運氣比較好,獲得了湖南省一等獎、全國二等獎的成績,就自己最大的感受而言,抛開團隊及其他因素,前期的積澱是比較重要的,數學模組化算法及應用這本書大部分内容要去弄懂,代碼盡量去實作,如果學校有教育訓練,盡量去參加,不要劃水,我所了解到的劃水的最後至多得了個省二省三,該付出時間的你一個都跑不掉。今年國賽我也報名了,希望我和各位都能取得一個滿意的成績。