永冻土层热传导问题数学建模

温度随时间，随距离变化的图像

图7 第二层温度随时间的变化图像

图8 第二层温度随距离的变化图像

17

图9 第二层温度随时间的变化图像

图10 第三层温度随距离的变化图像

18

图11 第四层温度随时间的变化图像

图12 第四层温度随距离的变化图像

问题二

模型二误差分析程序：

%% 求土壤温度分布 clear;clc; format short e

a=input(' 请输入系数a 的值：'); l=input(' 请输入长度l 的值：');

M=input(' 请输入将区间[0,l]等分的个数M：'); ot=input(' 请输入时间增量ot的值：'); n=input(' 请输入运行次数n 的值：'); ox=l/M;x0=zeros(M+1,1); for ii=1:M

x0(ii+1)=ii*ox; end

u=sin(pi*x0/l);%t=0 时u(x,t)的值 r=a^2*ot/(ox)^2; for ii=1:n %数据的输入

B=zeros(M-1,1);%存放系数矩阵主对角线元素

A=zeros (M-2,1);%存放系数矩阵主对角线元素下方次对角线的元素 C=zeros (M-2,1);%存放系数矩阵主对角线元素上方次对角线的元素 S=zeros(M-1,1);%存放右端的常数项 for ii=1:M-2

B(ii)=1+2*r;A(ii)=-r;C(ii)=-r; S(ii)=u(ii+1,1);

19

end

B(M-1)=1+2*r;S(M-1)=u(M,1);u(1,2)=0;u(M+1,2)=0; S(1,1)=S(1,1)+r*u(1,2);S(M-1,1)=S(M-1,1)+r*u(M+1,2); %追赶法

S(1)=S(1)/B(1);T=B(1);k=2; while k~=M

B(k-1)=C(k-1)/T;

T=B(k)-A(k-1)*B(k-1);

S(k)=(S(k)-A(k-1)*S(k-1))/T; k=k+1; end k=1;

while k~=M-1

S(M-1-k)=S(M-1-k)-B(M-1-k)*S(M-k); k=k+1; end

u(2:M,2)=S; %把结果放入矩阵u 中 u(:,1)=u(:,2);% 过河拆桥 end

%计算精确值，存放在u 的第二列 for x=0:M

u(x+1,2)=exp(-(pi*a/l)^2*n*ot)*sin(pi*x*ox/l); end

%计算最大相对误差 ez=zeros(M-1,1); for ii=2:M

ez(ii-1)=abs(u(ii,1)-u(ii,2))/u(ii,2); end

E=max(ez);

fprintf (' 最后时刻数值解与精确解分别为:\\n'); disp(u);

fprintf (' 差分法得到的结果与正确结果的最大相对误差为：'); disp([num2str(E*100) '%']); %画二维图比较

plot(x0,u(:,1),'g-',x0,u(:,2),'m*'); legend(' 数值解',' 精确解') xlabel('x'),ylabel('u(x,t)')

title(' 数值解与精确解比较') 问题三

温度含水量拟合程序 creatFit.m

20