一維熱傳導方程數值解法及matlab實現分離變量法和有限差分法

1、一維熱傳導方程的Matlab解法分離變量法和有限差分法問題描述有限長桿的熱傳導問題背景：一根長為 L 的均勻導熱細桿，側面絕熱，內部無熱源二其熱傳導系數為 k,比熱為 c,線密度為 d 求桿內溫度變化的規(guī)律。分析： 桿的溫度變化和熱里有關 設桿平躺在X軸上， 其端電在肝。 和x=L處。因為是細桿，且均勻導熱、制面絕荻，所以桿內熱里流動僅在K軸方向口因為均勻，所以k，C，。均為常數，考虛桿上AU到x+dx的一段（代表卜其而量為m=pdx,熱容量為cdm。設桿中的物流沿乂軸正向，強度為q（xt），溫度分布為u區(qū)t）。u0 0 xx+cxx+cLX物理定律，能:量守恒定律和熱傳導的Fourier定律

2、熱傳導的 Fourier 定律:若沿K方向有一定的溫度差， 在工方向也就有一定的熱量傳遞。從宏觀上看，單位時間內通過垂直Y方向的單位面積的熱量占溫度的空間變化率成正比.7duq=_#k（1）CXIFqfH熱流密度（強度），單位時間單位面積流過的熱量；上熱導率（熱傳導系數）能量守恒定律;因為內部無熱源，凈流入的熱量應該等于介質在此時間內溫度升高所需要的熱量。cdmdu=dQ=q(x.t)-q(x+dx.t)dt=q(x,t)dxdt/cdmdu=cpdxdu.cpdxdu=-qxdxdt=_%，即。外,=-9x(2)Q(X,Q(X,七)Q(x+dx,t)Q(x+dx,t)A2xx+dxxx+dx

3、介質內存在熱源時如果在介質內有熱量產生(:例如，有化學反應發(fā)生，或者通有電流，)，單位時間內單位體積介質產生的熱量吊尸&認因為熱傳導的Fourie定律沒有變化，所以仍然有,Su”質對于能量守恒定律，有cdmdu=-qT(x,t)dxdt+Fdxdt.即cWf+F(2)實驗原理分離變量法實驗原理有界長桿的熱傳導問題考察齊次熱傳導方程的混合問題（邊界條件都是第一類的情形）ut=(0 x0),、w(O,Z)=07z/(Z;Z)=0,u(xjy)(P(X),其中尹（約為給定的已知函數.下面用分離變量法（或稱駐波法）來求解定解問題（17）。首先令將其代入方程 2u,-QFK并分離變量得兩個常微分

4、方程右T”)=O.X”(H)+小工)=0：由邊界條件孤0/)=0,廿(/J)=0,可得X(0)=0,X(/)=0.(17)求邊值問題X”(x)+AY(x)=0,X(0)=X(l)=0.的非 0 解。(1)當 20 時，該問題有非平凡解。此時；M M、，(1)、4=4=()(刀=L2,.).居(x)=2sinW(=1；2,.).4現在考慮(。+加力，)=0,代入上方程得3+嚀尸 TG)=0,其通解為5=L2,).于是可得定解問題(17)中的一維熱傳導方程且且滿足齊將特征值(=L2.次邊界條件的具有變量分離形式的特解.nTDc/1Qane1sin,(18)n=ln=lI其中 4=B.g,是任意常數

5、。再利用初值條件(陽 0)=奴?？傻?.n7KI、&smsm丁=內)，AlI=o%=d*)sin與dx,(19)念、2k 式不)tnTDCu(x.t)=2anesinn n l/l/歡x)s山千及(18)(19)合在一起就是所求定解問題(17)的特解。%=a，。(0 x0),u(0,t)=0,u(J.t)=0?(1(丫 0)=o(x)7)若問題中的宮界條件出嬴第二類或第三類齊次邊界條件，解法類似。(18)(19)有限差分法、有限差分法的特點，有限差分方法（FDM）是計算機數值模擬最早采用的方法,至今仍被廣泛運用.該方法將求解域劃分為差分網格，用有限個網格節(jié)點代替連續(xù)的求解域.有限差分法

6、以Tayl口工級數展開等方法，把控制方程中的導數用網格節(jié)點上的函數值的差商代替送行離散，從而建立以網格節(jié)點上的值為未知數的代數方程組4該方法是一種直接將微分問題變?yōu)獒苋紗栴}的近似數值解法，數學概念直觀，表達簡單，是發(fā)展較早且比較成熟的數值方法。有限差分法的缺點是必需進行整個區(qū)域的剖分，并且要求網格比較規(guī)則，空間網格最好為直角網格。三,熱傳導方程(拋物方程)1.熱傳與方程的介紹為浮=優(yōu)Ta*0,0 x笈stdx1(x,0)=100,0 x0I*(1)用分離變量法和有限差分法解上面的問題。m(2)用分離變量法求出的u用。二也逼近，令m=100求出?1=1S由在一個坐標下畫出兩種方法求出的解和孔的圖

7、.(3)對上面求出的三個解，分別求出使得溫度低于50攝氏度的時間。(1)分離變量法(代碼)：x=0:0.1*pi:pi;y=0:0.04:1;x,t=meshgrid(x,y);s=0;m=length(j);%matlabfori=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=

8、pi構造一個1025列的矩陣(初始化為0)用于存放時間t和變量x橫坐標為x縱坐標為ts=(1/25)/(pi/10)A2;fprintf(穩(wěn)定性系數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);所得到的熱傳導圖形為:i分離

9、變量法(取前100項和)x=0:0.1*pi:pi;y=0:0.04:1;x,t=meshgrid(x,y);title(有限差分法解);s=0;fo門=1:100s=s+(200*(1-(-1)Ai)/(i*pi)*(sin(i*x).*exp(-iA2*t);end;surf(x,t,u);xlabel(x),ylabel(t),zlabel(T);title(分離變量法);axis(0pi010100);所得到的熱傳導圖形為：Ii有限差分法根據（1）我們有如下圖結論：比較可得這兩幅圖基本相同，有限差分法和分離變量法對本題都適應(3)第一種情況(取無窮項)在原來程序代碼白基礎上加上disp

10、(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);我們得到如下一組數據： clearclear127.324127.324122.3

11、315117.122.3315117.53485348112.9262103112.926210344983104.244044983104.2440100.100.1566156696.96.2294229492.456292.456230935.30935.3478347832.001332.00137S.7S.78597859當溫度低于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)

12、;r=0.04/(0.1*pi)A2;fprintf(穩(wěn)定性系數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);愚定性系數為：0.40530.405375.3030i75.3030i99.99.3634363472.72.45354535100.100.000069,6852000069,685299.99.932893

13、2867.003267.003299,731264.411299,731264.411298,98,9021902161.61.9101910197,399197,399159.59.4997499775,696775,696772.72.7286728669,69,8769876967H37067H37064.64.5045504561.61.9752975259.59.5451545157,210357,210354,54,9671967152,52,81188118StfTTStfTTJnbxJnbx碼至)46.46.8399839995,325095,325057.57.1737173

14、792.S37854.92.S37854.9450945090,082390,082352.52.79647964371719百 W 而、84,189384,1893kiki理)81J92646.S346.S3&6&678.78.22162216當溫度低于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)A2;fprintf(穩(wěn)定性系數S為:n);d

15、isp(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(有限差分法解);disp(u(5,:);得到如下結果Columns1through7Columns1through7100.0000100,0000100,0000100,000

16、097,3020100.0000100,0000100,0000100,000097,3020Colums8through14Colums8through1486.519682.59227E.7i5375.004686.519682.59227E.7i5375.0046我們知19歹U為50.350520列是數據為47.8902所以時間t為20*0.04=0.78結論：比較一二三種情況，我們得到不同的時間，這是由于:1、加和不同一種為100,一種為無窮；2、采用的方法不同：一種為分離變量法，一種為有限差分法造成的。第一題完94,164290.387894,164290.387871.404757

17、-952964.654671.404757-952964.6546ColumnsColumns15through2115through216E507858.509455.654152.936550.35056E507858.509455.654152.936550.350517,890245.549817,890245.5498Colwinis22through2522through2543.323541.205943.323541.205939.191937.276039.191937.2760作業(yè)3熱傳導方程的向前差分格式:H/T=sir-i+(1-2s)加+sid.Il-rl、zII-1

18、中保證穩(wěn)定性要使02o取恰當的s使得并用算出的u(i,j)數據說明取這個s時，數值方法不穩(wěn)定。解：根據課本知識：1,穩(wěn)定性判別Von-Neumann稔定性在判斷有限差分近似的穩(wěn)定性方法中，以Mon-N切nn方法使用較為廣泛，它僅適用于線性常系數的有限差分近似.其過程如下：首先，要研究的差分方程可寫為：X。渭黑=E-如，+1/Ic門u=su工+(1-2s)u+su.J7+1z7JTka2田二 K其次，對時進行變量分離*%=Zexp，平xj最后將;Vnexpzez代人所考察的有限差分方程。VnJrXexpzcrxy=sVnexpzofxy+14-(12s)Pexp7axexpfzcxL】=skex

19、pUc+(-2sVn+sVnexp-zAvn+1V7r5(cosahisina*=12s(1cosah). .2 2ah,11145sirr12從上式，我們看出，1,1Y- -2 2ah114sm214sin2ochT-下面說明，在什么條件下能使12,我們得到的數據就會不穩(wěn)定，所畫出的圖形失真也就會很大。 一下就S72,兩種情況進行討論(以第一題為例進行討論)當S12時，取步長k=1/25,h=pi/10,S=0.4053;u=zeros(10,25);%t=1x=pi10行25列橫坐標為x縱坐標為ts=(T25)/(pi/10)A2;fprintf(穩(wěn)定性系數S為:n);disp(s);fo

20、ri=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(有限差分法解s12,時取步長k=1/25,h=pi/20,S=1.6211u=zeros(20,25);%t=1x=pi20行25列橫坐標為x縱坐標為ts=(W25)/(pi/20)A2;fpri

21、ntf(穩(wěn)定性系數S為:n);disp(s)；fo門=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(有限差分法解s1/2時)；所得數據為Columns10through180 00 00 00 00 00 00 00 00 0-0.0000-0

22、.00000.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.-0.0000000 00.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.-0.00000000-0.-0.0000000 0-0.0000-0.00000.00000.0000-0.-0.000000000.00000.0000-0.-0.00000000-0.-0

23、.000000000.00000.0000-0.-0.000000000.00000.0000-0.0000-0.00000.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.-0.0000000 00.00000.0000-0.0000-0.00000.00000.0000-0.-0.000000000.00000.0000-0.0000-0.00000.00000.0000-0.0000-0.00000.00000.0000-0.-0.000000000.00

24、000.0000-0.0000-0.00000.00000.0000-0.0000-0.00000.00000.0000-0.0000-0.00000.00000.0000-0.000-0.0000 00.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.0000-0.00000.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.-0

25、.000000000.00000.0000-0.-0.0000000 00.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.0000-0.00000.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.-0.0000000 00.00000.0000-0.-0.000000000.0

26、0000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.0000-0.00000.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.-0.0000000 00.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.-

27、0.000000000.00000.0000-0.0000-0.00000.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.-0.0000000 00.00000.0000-0.0000-0.00000.00000.0000-0.0000-0.00000.00000.0000-0.0000-0.00000.00000.0000-0.0000-0.00000.00000.0000-0.0000-0.00000.00000.0000-0.-0.000000000.0

28、0000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.000-0.0000 00 00 00 00 00 00 00 00 00 00.00000.0000-0.-0.Q000Q0000.00000.0000-O.-O.00Q000Q00.0.00000000-O.OOOQ-O.OOOQ0.00000.0000_Q._Q.000000000.0.000000000.00000.0000Q.0000Q.0000-p.OQQQ-p.OQQQ0.0.00000000-0.-0.000000000i00000i0000-0-00

29、00-0-00000.00000.0000-0-0、000000000.0.000000000.0.000000000.00000.0000-a.-a.000000000.0.00000000-0.-0.000000000.0.ODDOODDO-0.-0.OQOOOQOO0.0.0000000 00.00000.0000o,o,oooooooo0.00000.00000.00000.0000rQ.rQ.000000000.00000.0000-o.oooo-o.oooo0.00000.0000T.QQQ0T.QQQ00,00000,00000.00000.00000.0.000000000.0

30、0000.0000-0.-0.000000000,00000,0000-0.0000-0.00000,00000,0000-0.0000-0.00000.0.000000000.0.000000000.00000.0000-O.-O.000000000.0.oooooooo-0.0000-0.00000.0. 口 口0000-0.-0.000000000.0.00000000O.OQQQO.OQQQ口,00000000-0.0000-0.00000,00000,0000-0,-0,00000000Q.QQQQQ.QQQQ-0.0000-0.00000,000Q0,000Q-0,0000-0,0

31、0000.00000.0000-0.-0.000000000,00000,0000-0.-0.000000000.0.00000000-0,0000-0,00000.00000.0000-0.0000-0.00000.00000.00000 00.0.00000000-0.-0.00000000O.O.00000000-0.-0.000000000.0.00000000-0.-0.000000000.00000.0000-0.-0.0000000 00.00000.0000-0.0000-0.00000.00000.0000-OH000000000,0,00000000-0.-0.004000

32、400.00000.0000-0.0000-0.00000,00000,00000,00000,00000.0.oooooooo-fl.oooo-fl.oooo0.0.000000001Q.1Q.000000000.00000.0000-0.-0.000000000.QOOO0.QOOO-0.-0.0000000 00.0.000000000.00000.00000.00000.0000-O.-O.000000000.0.00000000-0.-0.000000000.00000.0000-0.-0.000000000.0.000000000.00000.00000.00000.00000.0

33、0000.0000Q Q、000000007,00007,00000.00400.0040-0.-0.000000000.00000.0000-0.0000-0.00000.00000.00000.0.000000000.00000.0000(J.(J.000000000.0.OOQOOOQO-0.-0.000000000.00000.0000-0.-0.000000000.0.0000000 00.00000.0000o.ooooo.oooo0.00000.0000Q.Q.00000000-0.-0.00000000O.OOOQO.OOOQ-0.0000-0.00000.00000.0000

34、-0.-0.oooooooo0.00000.00000.00000.000000000000-0.-0.000000000.0.00000000-0.-0.000000000-00000-00007*00007*00000.0.000000000.0.000000000.0.00000000-0.-0.00000000O.O.00000000-0.-0.00000000O.O.00000000-0.-0.000000000.0.OQOOOQOO-0.-0.0000000 00.00000.0000-0,0000-0,00000.00000.0000-0.0000-0.00000.0.00000

35、000-O.OQQO-O.OQQOo.ooooo.oooo-0.0000-0.0000O.QQOOO.QQOO0 00 00 00 0Q Q0 00 00 00 0Columnathrough90.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.0000-0.00000.00000.0000-0.-0.000000000.00000.0000-0.-0.000000000.00000.0000-0.-0.000

36、000000.00000.0000-0.-0.0000000 00 00 00 00 00 00 00 00 00 0Columns19through25Columns19through250 00 00 00 0Q Q0 00 0O.OQOOO.OQOO0.00000.0000-0.0001-0.00010.00080.0008-0.0059-0.00590.03S40.03S4-0.2629-0.26290.00000.00000*00000*00000.00040.00040.0.002500250.01620.0162-0.-0.101510150.6190.6194 4-0.0000

37、-0.00000.00020.0002-0.0009-0.00090.00570.0057-0.0343-0.03430.20230.2023-L-L179E179E0,00010,0001-0,0003-0,00030.00190.0019-0,0107-0,01070.06110.0611-0.3465-0.3465I.9554I.9554-0.0001-0.00010.00060.0006-0.0031-0.00310,01710,0171-0.0949-0.09490,0,52475247-2.-2.896S896S0.00010.0001-0.OOQ8-0.OOQ80.0.00440

38、044-0.0241-0.02410.13130.1313-0.7147-0.71473.3.3903909 9-0.0002-0.0002040011040011-0.0057-0.00570.03050.03050.16440.16440*S8690*S869-4.7867-4.78670*00020*0002-0.0012-0.00120.0.00650065-0.0351-0.03510.18840.1884-L-L011301135.4345.4341 1-0,0002-0,00020,00130,0013-0,-0,006900690,03710,0371-0.1988-0.198

39、81.06641.0664-5.7255-5.72550.00020.0002-0.0013-0.00130,00670,0067-0.0362-0.03620.19450.1945-1,0455-1,04555,6255,6251 1-0.0002-0.00020.00110.0011-0i0061-0i0061Q.0328Q.0328-0.1771-0.17710.95730.9573-5.1756-5.17560.00020.0002-4.0009-4.00090.00510.0051-0.0277-0.02770.15110.1511-0.8230-0.82304,4,4804800 0-0.0001-0.0001。 ,00070007-0.0040-0.00400.02210.0221-0.1216-0.12160.66830.6683TSS93TSS930,00010,0001-0,0005-0,00050.00300.0030Y,0167Y,01670.09300.0930- -0.5100.5102,8562,8564 4-0.0001-0.00010.00040.0004-0.0022-0.00220