一维热传导方程数值解法及matlab实现分离变量法和有限差分法

问题描述有限长杆的热传导问题有限长杆的热传导问题背景：一根长为L的均匀导热细杆，侧面绝热，内部无热源n其热传导系数为k,比热为C,线密度为供求杆内温度变化的规律匚分析：杆的温度变化和热里有关。设杆平躺在其轴上，其端点在x=0和*=1_处。囱为是细杆，且均匀导热、例］面绝热，斫以杆内热里流动仅在X轴方向。因为均匀，所以k,c,P均为常数。考虑杆上从品到乂+七的一段（代表），其质里为巾=pdx-热容里为匚dm。设杆中的热流沿*轴正向，强度为q-t），温度分布为u（xFth物理定律；・能量守恒定律和热传导的F口而立定律热传导的Fourier定律:若沿工方向有一定的温度差,在/方向也就有一定的热量传递。从宏观上看，单位时间内通过垂直t方向的单位面积的热量4与温度的空间变化率成正比。热流密度（强度），单位时间单位面积流过的热量美展热导率（热传导系数）

能量守恒定律二因为内部无热源，净流入的热量应该等于介质在此时间内温度升高所需要的热量。cdmdu=dQ=[q(x.t)-q(x+dx.t)]dt="qs(x?t)dxdtcdmdu=cpdxdit.cpdxdu=-q^dxdt

介质内存在热源时•如果在血质内有热量产生（例如，有化学反应发生，或者通有电流，•£）,单位时间内单•位体和介质产生的热量为尸a3因为热传导的Fourie定律没有变化，所以仍然有对于能量守恒定律，有cdmdu=—qx(x.t)dxdt-Fdxdt;Hpcit=-(2)实验原理分离变量法实验原理有界长杆的热传导问题一、考察齐次热传导方程的混合问题（边界条件都是第一类的情形）%= （0<X</.t>。），W仪（。］）=。：N&t）=0?"（阳0）=口（工），其中河行为给定的已知函数.下面用分离变量法（或称驻波法）来求解定解问题（孙首先令将其代入方程 工IL=U?%并分离变量得两个常微分方程「«）+勿叼（。=&X”（x）+M（x）=0,由边界条件以o,0=由祖界f）二。,可得阳。）=5x（r）=o.

求边值问题zu(x)+租x)=0, X(0)=X(l)=0.的非。解.(D当a<o时，该问题没有非平凡解:(2)当彳=0时，该问题也没有非平凡解二(3)当2>0时，该问题有非平凡解。此时鼻=&=(等尸(耳=12…).xQ)=5.sm牛5=上2t现在考虑 「(0+加之丁。)=0,将特征值后z=(尹5=124低入上方程得其通解为骞⑺工力，5=1二…)一于是可得定解问题(17)中的一维热传导方程且且满足齐次边界条件的具有变量分离形式的特解

竿，(18)TOC\o"1-5"\h\zw=l t4=:(飙>)sin詈双 (19)(18)(19)合在一起就是所求定解问题(17)的特解。巩=a2iiyx(0<x<Lt>0),

t笈(0乃二。？〃(『")=£ (1Iz/(x30)=火幻， ,)若问题中的边界条件出现第二类或第三类齐次边界条件，解法类似.有限差分法'有限差分法的特点・有限差分方法审明)是计算机数值模拟最早采用的方法；至今狷被广修运用。该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域中有限差分法以T的lu级数展开等方法:；把控制方程中的导数用网格节点上的函数值的差商代替进行离散「以而建立以网格节点上的值为未知数的代数方程组。・该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观;1•表达简单，是发展较早且比较成熟的数值方法；.一・有限差分法的缺点是面需进行整个区域的剖分：并且要求网格比较规则，空间网格最好为直角网格。J热传导方程(抛物方程)1热传导方程的分绍 [■加,居豆一/短«虱0j)=试上5)=。

«fc0)=/(.x)2.离散化ui=u(O.jij=0 u;'-=u(Lj)=0uf=u(M=f(jh)=ft(1)向前差分格式：k *•可以证明:当0"工!时，上述差分格式是稳定的口所以X的步长h和t的步长k取法要恰当口

实验目的利用分离变量法和有限差分法解热传导方程问题利

用matlab进行建模构建图形研究不同的情况下采用

何种方法从更深层次上理解热量分布与时间、空间

分布关系。模拟与仿真作业du_d2u

dtdx2<〃(苞0)=100,0<=〃(1/)=0,/>0(1)用分离变量法和有限差分法解上面的问题.m(2)用分离变量法求出的u用%二£4 逼近，令m=10。求出?2=1为■在一个坐标下画出两种方法求出的解和S刑的图.(3)对上面求出的三个解，分别求出使得温度低于50摄氏度的时间1(1)分离变量法(代码)：x=0:0.1*pi:pi;y=0:0.04:1;[x,t]=meshgrid(x,y);s=0;m=length(j);%matlab可计算的最大数相当于无穷fori=1:ms=s+(200*(1-(-1)Ai))/(i*pi)*(sin(i*x).*exp(-iA2*t));end;surf(x,t,s);xlabel('x'),ylabel('t'),zlabel('T');title('分离变量法(无穷)');axis([0pi010100]);所得到的三维热传导图形为：

有限差分法：u=zeros(10,25); %t=1x=pi构造一个1025列的矩阵(初始化为0)用于存放时间t和变量x横坐标为x纵坐标为ts=(1/25)/(pi/10)八2;fprintf('稳定性系数S为：\n');disp(s);fori=2：9u(i,1)=100;end;forj=1：25u(1,j)=0;u(10,j)=0;end;forj=1：24fori=2：9u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j);endenddisp(u);[x,t]=meshgrid(1:25,1:10);surf(x,t,u);xlabel('t'),ylabel('x'),zlabel('T');title('有限差分法解,)；所得到的热传导图形为:i分离变量法（取前100项和）x=0:0.1*pi:pi;y=0:0.04:1;[x,t]=meshgrid(x,y);s=0;fori=1:100s=s+(200*(1-(-1)入i))/(i*pi)*(sin(i*x).*exp(-i八2*t));end;surf(x,t,u);xlabel('x'),ylabel('t'),zlabel('T');title('分离变量法')；axis([0pi010100]);所得到的热传导图形为：Ii有限差分法根据（1）我们有如下图

(3)第一种情况(取无穷项)：在原来程序代码的基础上加上disp(s(:,6));可得出第六列(即x=pi/2)处温度随时间的变化情况x=0:0.1*pi:pi;y=0:0.04:1;[x,t]=meshgrid(x,y);s=0;m=length(j);%matlab可计算的最大数，相当于无穷fori=1:ms=s+(200*(1-(-1)Ai))/(i*pi)*(sin(i*x).*exp(-iA2*t));end;surf(x,t,s);

xlabel('x'),ylabel('t'),zlabel('T');title('分离变量法(无穷)’);axis([0pi010100]);disp(s(:,6));我们得到如下一组数据：75.696772.728669.75.696772.728669.876967.137064.504561.975259.545157.210354.967152.811846.8399127.324|0122.3315117.5348112.9262108.4983104.2440100.156696.229492.456288.830985.347882.001378.7859当温度低于50度是时间为t=23.5*0.04=0.94第二种情况（取前100项和）在原来程序代码的基础上加上disp(s(:,6));可得出第六列(即x=pi/2)处温度随时间的变化情况x=0:0.1*pi:pi;y=0:0.04:1;[x,t]=meshgrid(x,y);r=0.04/(0.1*pi)A2;fprintf('稳定性系数S为：’)disp(r);s=0;fori=1:100s=s+(200*(1-(-1)Ai))/(i*pi)*(sin(i*x).*exp(-iA2*t));end;surf(x,t,s);xlabel('x'),ylabel('t'),zlabel('T');title('分离变量法');axis([0pi010100]);disp(s(:,6));息定性系数s为：0.405375.303099.363472.4539100.000069.685299.982867.003299.731264.411298.902161.910197.399159.499795.325057.178792.837854.945090.082352.796487.1719&<730384.189381.192646.834678.2216当温度低于50度是时间为t=23.5*0.04=0.94第三种情况（有限差分法）在原来程序代码的基础上加上disp(u(5,:));可得出第五行(即x=pi/2)处温度随时间的变化情况u=zeros(10,25); %t=1x=pi10行25歹列横坐标为x纵坐标为ts=(1/25)/(pi/10)/2;fprintf('稳定性系数S为：\n');disp(s);fori=2：9u(i,1)=100;end;forj=1：25u(1,j)=0;u(10,j)=0;end;forj=1：24fori=2：9u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j);end

enddisp(u);[x,t]=meshgrid(1:25,1:10);surf(x,t,u);xlabel('t'),ylabel('x'),zlabel('T');title('有限差分法解,)；disp(u(5,:));得到如下结果Columns1through7100.0000100.0000100.0000100.000097.302094.1642Columns1through7100.0000100.0000100.0000100.000097.302094.164290.3878Columns8through14Columns8through1486.519682.592278.745375.004671.404767.952986.519682.592278.745375.004671.404767.952964.6546Columns15through21Columns15through2161.507858.509455.654152.936550.350547.890245.549861.507858.509455.654152.936550.350547.890245.5498Columns22through25Columns22through2543.323541.2059 39.1918 37.2760我们知19列为50.3505 20列是数据为47.8902所以时间t为20*0.04=0.78结论：比较一二三种情况，我们得到不同的时间，这是由于：1、加和不同一种为100，一种为无穷；2、采用的方法不同：一种为分离变量法，一种为有限差分法造成的。第一题完作业3热传导方程的向前差分格式:U/=S匕1+(1-2sM+SUi-X中保证稳定性要使。取恰当的S使得并用算出的u(i,j)数据说明取这个S时，数值方法不稳定解：根据课本知识：1.稳定性判别Von-Neumann稳定性在判断有限差分近似的稳定性方法中，以Von-Neumann方法使用较为广泛，它仅适用于线性常系数的有限差分近似口其过程如下:首先，要研究的差分方程可写为：+(l—2s)〃；+sUj_r其次,对明进行变量分离:最后将〃；HVKexp[iaXj]代人所考察的有限差分方程0J7向exppax/=sV"exp[zcrxy+1]+(1-2s)V"exp[zYzx7]+sV/7expfztzx^]定义为

放大系

数Vn+Y=s，"cxp[zU川+(1-2sWn+sVnexp[-iah]vn+1=l-45sin2-^-下面说明，在什么条件下能使<1对所有的a成立「-l<l-45sin2—<1

2从上式，我们看出,-1<1-45sin2—21我们知道稳定性系数S是衡量由差分格式所得的数据是否稳定的因数，如果S>12,我们得到的数据就会不稳定，所画出的图形失真也就会很大。一下就S<1/2和S>1/2,两种情况进行讨论(以第一题为例进行讨论)当S<S时，取步长k=1/25,h=pi/10,S=0.4053;u=zeros(10,25);%t=1x=pi10行25歹列横坐标为x纵坐标为ts=(1/25)/(pi/10)八2;fprintf('稳定性系数S为:\n');disp(s);fori=2:9u(i,1)=100;end;forj=1:25u(1,j)=0;u(10,j)=0;end;forj=1:24fori=2:9u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j);endenddisp(u);[x,t]=meshgrid(1:25,1:10);surf(x,t,u);xlabel('t'),ylabel('x'),zlabel('T');title('有限差分法解s<1/2时')；获得的数据为

稳定性系数5为:0.4053Columns1through9000000000100.000059.471551.794243.632339.624235.951833.418131.212529.3941100.0000100.000083.574477.351470.187265.651961.394157.938454.7629100.0000100.0000100.000093.343089.559884.846280.843576.829173.1004100.0000100.0000100.0000100.000097.302094.164290.387886.519682.5922100.0000100.0000100.0000100.000097.302094.164290.387886.519682.5922100.0000100.0000100.000093.343089.559834.846280.843576.829173.1004100.0000100.000083.574477.351470.187265.651961.394157.938454.7629100.000059.471551.794243.632339.624235.951833.418131.212529.3941000000000CoIlbtitls10through1800000000027.762726.298724.945823.638022.506521.392620.33S619.339518.391251.913249.259346.788044.460742.264540.184538.212036.339134.559869.515466.122262.S87859.814656.890254.109451.464348.943646.555878.745375.004671.404767.952964.654661.507858.509455.654152.936578.745375.004671.404767.952964.654661.507858.509455.654152.936569.515466.122262.S87859.814656.890254.109451.464348.943646.555851.913249.259346.788044.460742.264540.184538.212036.339134.559827.762726.298724.945823.638022.506521.392620.338619.339518.3912000000000Colmms19through25000000017.490416.634415.820715.046914.311213.611512.946032.868731.261029.732228.278426.895925.581024.330444.280042.115540.056738.098636.236234.464932.780150.350547.890245.549843.323541.205939.191837.276050.350547.890245.549843.323541.205939.191837.276044.280042.115540.056738.098636.236234.464932.780132.868731.261029.732228.278426.895925.581024.330417.490416.634415.820715.046914.311213.611512.94600000000所得到的图形为

fprintf(稳定性系数S为:\n');disp(s);fori=2:19u(i,1)=100;end;forj=1:25u(1,j)=0;u(10,j)=0;end;forj=1:24fori=2:19u(i,j+1)=s*u(i+1,j)+(1-2*s)*u(i,j)+s*u(i-1,j);endenddisp(u);[x,t]=meshgrid(1:25,1:20);surf(x,t,u);xlabel('t'),ylabel('x'),zlabel('T');title('有限差分法解s>1/2时')；

所得数据为Colunins1through90000000000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.00000.00000.0000-|o.00000.0000-0.00000.0000-0.00000.0000-0.00000.00000.00000.0000-0.00000.0000-0.00000.0000-0.00000.00000.00000.00000.00000.0000-0.00000.0000-0.00000.0000-0.00000.00000.00000.00000.0000-0.00000.0000-0.00000.0000-0.00000.00000.00000.0000-0.00000.0000-0.00000.0000-0.00000.00000.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.000000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.00000.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.00000.00000.0000-0.00000.0000-0.00000.0000-0.00000.00000.00000.00000.00000.0000-0.00000.0000-0.00000.0000-0.00000.00000.00000.00000.00000.0000-0.00000.0000-0.00000.00000.00000.00000.00000.0000-0.00000.0000-0.00000.0000-0.00000.00000.00000.0000-0.00000.0000-0.00000.0000-0.00000.00000.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000000000000CoIujtliis10throug0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.0000-0.0000-0.00000.0000-0.00000.0000-0.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.00000.0000-0.0000000000000

Columns19through250Q00Q000.00000.0000-0.00010.0008-0.00590.0394-0.25290.0000-0.00000.0004-0.00250.0162-0.10150.6194-0.00000.0002-0.00090.0057-0.03430.2023-1.17980.0001-0.00030.0019-0.01070.0611-0.34651.9554-0.00010.0006-0.00310.0171-0.09490.5247-2.89680.0001-0.00080.0044-0.02410.1313-0.71473.8909-0.00020.0011-0.00570.0305-0.16440.8869-4.78670.0002-0.00120.0066-0.03510.1884-1.01135.4341-0.00020.0013-0.00690.0371-0.19881.0664-5.72550.0002-0.00130.0067-0.03620.1945-1.04555.6251-0.00020.0011-0.00610.0328-0.17710.9573-5.17560.0002-0.00090.0051-0.02/70.1511-0.82304.4808-0.00010.0007-0.00400.0221-0.12160.6683-3.66930.0001-0.0