偏微分一维热传导问题

偏微分大作业一维热传导方程问题运用隐式格式求解数值解目录问题描述31 解析解分离变量法32 数值解隐式格式53 证明隐式格式的相容性与稳定性54 数值解分析与Matlab实现65 数值解与解析解的比较96 随时间变化的细杆上的温度分布情况117稳定后细杆上的温度分布情况12参考文献13附录14有限长杆的一维热传导问题问题描述一根单位长度的细杆放入100的沸水中，当细杆的温度达到100时取出。假设细杆四周绝热；在时间t=0时，细杆两端浸入0的冰水中。一维热传导方程：，现在令，从而可知本题：。现在要求细杆温度分布：。1 解析解分离变量法热传导偏微分方程： (1) 其中， 首先令： (2)将(2)式带入(1)式得： 于是可得： 可以得到两个微分方程：先求解空间项：当时， 由于可知:由于解的收敛性， 则此时是平庸解。当时， 则此时是平庸解。当时， ，其中。 所以， 因为所以， 则，初始条件： 最终， 2 数值解隐式格式目前，研究热传导问题特别是非稳态热传导问题十分重要。这里使用隐式格式。利用，关于t进行向前差商： ；关于x进行二阶中心差：；代入偏微分方程可以得到隐式差分格式： （1）3 证明隐式格式的相容性与稳定性（1）相容性 代入隐式格式得： (2) 将（2）与原微分方程相减，得到截断误差： 所以此隐式格式与原微分方程相容。（2）稳定性 令网格比为，则可以将（1）式改写得到： （3）首先令： （4）将（4）代入（3）式，根据欧拉公式化简得： （5）故得放大因子是：所以根据Fourier方法，隐式格式恒稳定。4 数值解分析与Matlab实现（1） 边值与初值离散化将边值与初值离散化，与式（3）联立得差分线性方程组： ， ， 再将方程组改写成的形式：本题的边界条件均为零。所以可以将上式改写。（2） Matlab的实现 杆长1米，时间2秒。设计空间步长h=0.1和时间步长t=0.01，网格比是。从而得到划分的空间网格点数是M1+1，时间网格点数是M2+1。先设初始的温度矩阵U(M2+1,M1+1)。再将边界条件和初始条件编写到表示温度分布的矩阵中。具体代码可见最后附录。 编写矩阵A核心代码：对角线：A(i,i) = 1+2r 对角线的右方和下方：A(i,i+1) = -r; A(i+1,i) = -r; 下面就要运用进行迭代。当k=1时，A*U(2,j)=U(1,j)当k=2时，A*U(3,j)=U(2,j)当k=3时，A*U(4,j)=U(3,j)以此迭代下去直到k=M2。就可以得到整个温度随时间和空间的分布矩阵U。 数值解画图，如图1(a)和图1(b)所示。图1(a) 数值解的温度分布图现在将着色平稳过渡。图1(b) 着色平稳过渡的数值解的温度分布图5 数值解与解析解的比较 首先，我们需要将解析解离散化，解析解中有一项，当n越来越大时，会快速趋于0，故我们可以取n=8000。现在来证明可行性，在matlab里的工作空间运算。 将解析解的温度分布画出来，数值解画图，如图2所示。图2 解析解的温度分布图将数值解与解析解相减，得到误差图。如图3(a)和图3(b)，我们从图3(a)上可以看出空间上的误差，在边界处误差比较大。图3（a） 数值解与解析解空间误差我们从图3(a)上可以看出时间的误差，在时间的最开始，处误差最大，然后又有一个小的波动，最后就误差渐渐变小，最后趋于0。图3（b） 数值解与解析解时间误差6 随时间变化的细杆上的温度分布情况从数值解的温度分布三维图，如图4(a)和图4(b)可以看出随着时间的增加，细杆温度下降最后趋于0。从物理角度来说：细杆的温度会不断地向两端扩散，热量会慢慢散失，最终随着时间的增加，细杆的温度会趋于0。图4(a) 细杆温度随时间的变化图现取细杆中心处一点，观看它随时间的温度变化情况。图4(b) 细杆中央（x=0.5）温度随时间的变化图7稳定后细杆上的温度分布情况从图像上可以看出，最后稳定的情况下，细杆的温度是0。参考文献1 冯立伟热传导方程几种差分格式的MATLAB数值解法的比较J沈阳化工大学，辽宁沈阳2011(6)2 一维热传导方程数值解法及Matlab实现EB/OL2014-11-20http:/www.doc88.com/p-65221891945.html附录代码：%此程序用于解决一维热传导方程:ut-a2uxx = 0 %边界条件：u(0,t) = u(L,t) = 0 % %初始条件：u(x,0) = 100, x!=0和L % u(0,0) = 0 % u(L,0) = 0 %其中,a2 = 1, L = 1 %clc;clear all;%区域及划分网格L = 1; %单位长度的细杆T = 2; %时间h = 0.1; % 空间的划分 %t = 0.01; % 时间的划分%r = t/(h*h); %网格比 %设计步长M1 = L/h;M2 = T/t; %构造边界条件 %构造的矩阵：U(时间，空间）U = zeros(M2+1,M1+1); %编程包含边值，如U(k,1)=u(0,t)for k = 1:M2+1 %时间划分了M2份，有M2+1个节点 U(k,1) = 0; %两个边界处温度恒为零 U(k,M1+1) = 0; end; %构造初始条件for j = 2:M1 %位置划分了M1份，有M1+1个节点 U(1,j) = 100;end; U(1,1) = 0; U(1,M1+1) = 0; %差分格式的矩阵形式 A*U(k+1,j)=U(k,j)%构造矩阵AA = zeros(M1-1);for i = 1:M1-1 A(i,i) = 1+2*r;end;for i = 1:M1-2 A(i,i+1) = -r; A(i+1,i) = -r;end; %构造AU=B中的B %本题边值的特殊，矩阵B大大简化了B = zeros(M1-1,1);for k = 1:M2 j = 2:M1; B(j-1,1) = U(k,j); x = AB; for j = 2:M1 U(k+1,j) = x(j-1); %k+1时刻的不同位置的温度分布 end;end; %作图x = 0:h:1;y = 0:t:2;xx,yy=meshgrid(x,y);figure(1);surf(xx,yy,U);shading flattitle(一维热传导方程-数值解-温度分布图);xlabel(位置x);ylabel(时间t);zlabel(温度T); figure(2)s = 0;for i= 1:8000 s = s+(200*(1-(-1)i)/(i*pi)*sin(i*pi*xx).*exp(-i2*pi2*yy);end;surf(xx,yy,s);title(一维热传导方程-解析解-温度分布图);xlabel(位置x);ylabel(时间t);zlabel(温度T); figure(3)x = 0:h:1;y = 0:t:2;xx,yy = meshgrid(x,y);dd = U-s;surf(xx,yy,dd);title(一维热传导方程-误差-温度分布图);xlabel(位置x);ylabel(时间t);zlabel(误差（数值解减解析解）); figure(4)z = zeros(M2+1,1);if mod(M1+1),2)=0 i = 1:M2+1; z(i,1) = U(i,M1/2);else i = 1:M2+1; z(i,1) = U(i,(M1+1)/2);end;plot(z);title(温度随时间增加的趋势图);xlabel(时间t);ylabel(温度T);17