1、偏微分大作业一维热传导方程问题运用隐式格式求解数值解目录问题描述 .31 解析解 分离变量法 .32 数值解 隐式格式 .53 证明隐式格式的相容性与稳定性 .54 数值解 分析与 Matlab 实现 .65 数值解与解析解的比较 .96 随时间变化的细杆上的温度分布情况 .117 稳定后细杆上的温度分布情况 .12参考文献 .13附录 .14有限长杆的一维热传导问题问题描述一根单位长度的细杆放入 100的沸水中，当细杆的温度达到 100时取出。假设细杆四周绝热；在时间 t=0 时，细杆两端浸入 0的冰水中。一维热传导方程： ，现在令 ，从而可知本题： 。现20txua21a0txu在要求细杆
2、温度分布： 。(,)ut1 解析解分离变量法热传导偏微分方程：0txu(1)(,)(1,t)，其中，0x， 或()x1(,)，首先令： (2)(,)uxtXxTt将(2)式带入(1)式得：(0于是可得：)(txTX可以得到两个微分方程： )0tt()x先求解空间项：当 时， 0)xXAeB4由于 (0,t)1,t0.ut可知:由于解的收敛性， B=0XAeA则此时是平庸解。当 时， 0()x1,B则此时是平庸解。当 时， ，其中 。()cosinXxAkxk00(1)i,12,3.B所以， ， sin()xx因为 2n所以， ， 2()ntTtCe1,23.则， 21,si()ntuxtDx初
3、始条件： (0)(，1sin)(nxx， 02()nDd1si)(co(1cosxnn200lim=)nD当 时 ，最终，5， 2120(,)()sin()ntnuxtex1,23.2 数值解隐式格式目前，研究热传导问题特别是非稳态热传导问题十分重要。这里使用隐式格式 。1利用 ，关于 t 进行向前差商： ；关于 x 进行二阶中心(,)uxt 1kjjUtA差： ；112()kkjjjUA代入偏微分方程可以得到隐式差分格式：1112()kkkjjjjjtxA（1）3 证明隐式格式的相容性与稳定性（1）相容性 21 223312 2Taylor1=+()1()1kjjkkjjkkjjUtttUx
4、txtt AAAAA根 据 展 开 ：代入隐式格式得：(2)222()=+()UUttxt xA将（2）与原微分方程相减，得到截断误差 ：621=-()0UtxA所以此隐式格式与原微分方程相容。（2）稳定性令网格比为 ，则可以将（1）式改写得到： 2rtxA11()kkkkj jjjUrU（3）首先令：1(-11(+1+kkIjjjjkIjjkjjeUe） ）（4）将（4）代入（3）式，根据欧拉公式化简得：12cos)k krU（5）故得放大因子是： 11+2(cos)kGr所以根据 Fourier 方法，隐式格式恒稳定。4 数值解分析与 Matlab 实现（1） 边值与初值离散化将边值与初值
5、离散化，与式（3）联立得差分线性方程组：， 111(2)kkkkj jjjrUrU(0,12,M-)jN70=(),jjUx(0,12,M)jkNk， L(,)再将方程组改写成 的形式：AB110223312211(1)M1+2 +122kkkkkkMr Urrrrr 本题的边界条件均为零。所以可以将上式改写。 1122331221(1)M1+2122kkkkMr Urr （2） Matlab 的实现 杆长 1 米，时间 2 秒。设计空间步长 h=0.1 和时间步长 t=0.01，网格比是 。2trh从而得到划分的空间网格点数是 M1+1，时间网格点数是 M2+1。先设初始的温度矩阵 U(M2
6、+1,M1+1)。再将边界条件和初始条件编写到表示温度分布的矩阵中。具体代码可见最后附录。 编写矩阵 A核心代码：对角线：A(i,i) = 1+2r对角线的右方和下方：A(i,i+1) = -r;A(i+1,i) = -r;8 下面就要运用 进行迭代。*(1,)(,)AUkjkj当 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) 数值解的温度分布图现在将着色平稳过渡。9图 1(b) 着色平稳过渡的数值解的温度分布图5 数值解与解析解的比较 首先，我们需要将解析解离散化，解析解中有一项 ，当 n 越来越2nte大时，会快速趋于 0，故我们可以取 n=8000。现在来证明可行性，在matlab 里的工作空间运算。将解析解的温度分布画出来，数值解画图 ，如图 2 所示。210图 2 解析解的温度分布图将数值解与解析解相减，得到误差图。如图 3(a)和图 3(b)，我们从图 3(a)上可以看出空间上的误差，在边界处误差比较大。
