一維熱傳導(dǎo)方程數(shù)值解法及matlab實(shí)現(xiàn)分離變量法和有限差分法

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

2、熱傳導(dǎo)的 Fourier 定律:若沿K方向有一定的溫度差， 在工方向也就有一定的熱量傳遞。從宏觀上看，單位時(shí)間內(nèi)通過垂直Y方向的單位面積的熱量占溫度的空間變化率成正比.7duq=_#k（1）CXIFqfH熱流密度（強(qiáng)度），單位時(shí)間單位面積流過的熱量；上熱導(dǎo)率（熱傳導(dǎo)系數(shù)）能量守恒定律;因?yàn)閮?nèi)部無熱源，凈流入的熱量應(yīng)該等于介質(zhì)在此時(shí)間內(nèi)溫度升高所需要的熱量。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、介質(zhì)內(nèi)存在熱源時(shí)如果在介質(zhì)內(nèi)有熱量產(chǎn)生(:例如，有化學(xué)反應(yīng)發(fā)生，或者通有電流，)，單位時(shí)間內(nèi)單位體積介質(zhì)產(chǎn)生的熱量吊尸&認(rèn)因?yàn)闊醾鲗?dǎo)的Fourie定律沒有變化，所以仍然有,Su”質(zhì)對于能量守恒定律，有cdmdu=-qT(x,t)dxdt+Fdxdt.即cWf+F(2)實(shí)驗(yàn)原理分離變量法實(shí)驗(yàn)原理有界長桿的熱傳導(dǎo)問題考察齊次熱傳導(dǎo)方程的混合問題（邊界條件都是第一類的情形）ut=(0 x0),、w(O,Z)=07z/(Z;Z)=0,u(xjy)(P(X),其中尹（約為給定的已知函數(shù).下面用分離變量法（或稱駐波法）來求解定解問題（17）。首先令將其代入方程 2u,-QFK并分離變量得兩個(gè)常微分

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)當(dāng) 20 時(shí)，該問題有非平凡解。此時(shí)；M M、，(1)、4=4=()(刀=L2,.).居(x)=2sinW(=1；2,.).4現(xiàn)在考慮(。+加力，)=0,代入上方程得3+嚀尸 TG)=0,其通解為5=L2,).于是可得定解問題(17)中的一維熱傳導(dǎo)方程且且滿足齊將特征值(=L2.次邊界條件的具有變量分離形式的特解.nTDc/1Qane1sin,(18)n=ln=lI其中 4=B.g,是任意常數(shù)

5、。再利用初值條件(陽 0)=奴。可得2.n7KI、&smsm丁=內(nèi))，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)有限差分法、有限差分法的特點(diǎn)，有限差分方法（FDM）是計(jì)算機(jī)數(shù)值模擬最早采用的方法,至今仍被廣泛運(yùn)用.該方法將求解域劃分為差分網(wǎng)格，用有限個(gè)網(wǎng)格節(jié)點(diǎn)代替連續(xù)的求解域.有限差分法

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

7、.(3)對上面求出的三個(gè)解，分別求出使得溫度低于50攝氏度的時(shí)間。(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);所得到的三維熱傳導(dǎo)圖形為：可計(jì)算的最大數(shù)相當(dāng)于無窮有限差分法：u=zeros(10,25);%t=1x=

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

9、變量法(取前100項(xiàng)和)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);所得到的熱傳導(dǎo)圖形為：Ii有限差分法根據(jù)（1）我們有如下圖結(jié)論：比較可得這兩幅圖基本相同，有限差分法和分離變量法對本題都適應(yīng)(3)第一種情況(取無窮項(xiàng))在原來程序代碼白基礎(chǔ)上加上disp

10、(s(:,6);可得出第六列(即x=pi/2)處溫度隨時(shí)間的變化x=0:0.1*pi:pi;y=0:0.04:1;x,t=meshgrid(x,y);s=0;m=length(j);%matlab可計(jì)算的最大數(shù)，相當(dāng)于無窮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);我們得到如下一組數(shù)據(jù)： 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當(dāng)溫度低于50度是時(shí)間為t=23.5*0.04=0.94第二種情況(取前100項(xiàng)和)在原來程序代碼的基礎(chǔ)上加上disp(s(:,6);可得出第六列(即x=pi/2)處溫度隨時(shí)間的變化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)定性系數(shù)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);愚定性系數(shù)為：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當(dāng)溫度低于50度是時(shí)間為t=23.5*0.04=0.94第三種情況（有限差分法）在原來程序代碼的基礎(chǔ)上加上disp(u(5,:);可得出第五行(即x=pi/2)處溫度隨時(shí)間的變化情況u=zeros(10,25);%t=1x=pi10行25列橫坐標(biāo)為x縱坐標(biāo)為ts=(1/25)/(pi/10)A2;fprintf(穩(wěn)定性系數(shù)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,:);得到如下結(jié)果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列是數(shù)據(jù)為47.8902所以時(shí)間t為20*0.04=0.78結(jié)論：比較一二三種情況，我們得到不同的時(shí)間，這是由于: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熱傳導(dǎo)方程的向前差分格式:H/T=sir-i+(1-2s)加+sid.Il-rl、zII-1

18、中保證穩(wěn)定性要使02o取恰當(dāng)?shù)膕使得并用算出的u(i,j)數(shù)據(jù)說明取這個(gè)s時(shí)，數(shù)值方法不穩(wěn)定。解：根據(jù)課本知識：1,穩(wěn)定性判別Von-Neumann稔定性在判斷有限差分近似的穩(wěn)定性方法中，以Mon-N切nn方法使用較為廣泛，它僅適用于線性常系數(shù)的有限差分近似.其過程如下：首先，要研究的差分方程可寫為：X。渭黑=E-如，+1/Ic門u=su工+(1-2s)u+su.J7+1z7JTka2田二 K其次，對時(shí)進(jìn)行變量分離*%=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,我們得到的數(shù)據(jù)就會(huì)不穩(wěn)定，所畫出的圖形失真也就會(huì)很大。 一下就S72,兩種情況進(jìn)行討論(以第一題為例進(jìn)行討論)當(dāng)S12時(shí)，取步長k=1/25,h=pi/10,S=0.4053;u=zeros(10,25);%t=1x=pi10行25列橫坐標(biāo)為x縱坐標(biāo)為ts=(T25)/(pi/10)A2;fprintf(穩(wěn)定性系數(shù)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,時(shí)取步長k=1/25,h=pi/20,S=1.6211u=zeros(20,25);%t=1x=pi20行25列橫坐標(biāo)為x縱坐標(biāo)為ts=(W25)/(pi/20)A2;fpri

21、ntf(穩(wěn)定性系數(shù)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時(shí))；所得數(shù)據(jù)為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