偏微分一維熱傳導(dǎo)問題

偏微分一維熱傳導(dǎo)問題偏微分一維熱傳導(dǎo)問題偏微分一維熱傳導(dǎo)問題偏微分一維熱傳導(dǎo)問題編制僅供參考審核批準(zhǔn)生效日期地址：電話：傳真:郵編:偏微分大作業(yè)一維熱傳導(dǎo)方程問題——運(yùn)用隱式格式求解數(shù)值解目錄問題描述 31解析解——分離變量法 32數(shù)值解——隱式格式 53證明隱式格式的相容性與穩(wěn)定性 54數(shù)值解——分析與Matlab實(shí)現(xiàn) 65數(shù)值解與解析解的比較 96隨時(shí)間變化的細(xì)桿上的溫度分布情況 117穩(wěn)定后細(xì)桿上的溫度分布情況 12參考文獻(xiàn) 13附錄 14有限長桿的一維熱傳導(dǎo)問題問題描述一根單位長度的細(xì)桿放入100℃的沸水中，當(dāng)細(xì)桿的溫度達(dá)到100℃時(shí)取出。假設(shè)細(xì)桿四周絕熱；在時(shí)間t=0時(shí)，細(xì)桿兩端浸入0℃的冰水中。一維熱傳導(dǎo)方程：，現(xiàn)在令，從而可知本題：?，F(xiàn)在要求細(xì)桿溫度分布：。1解析解——分離變量法熱傳導(dǎo)偏微分方程：(1)其中，首先令：(2)將(2)式帶入(1)式得：于是可得：可以得到兩個(gè)微分方程：先求解空間項(xiàng)：當(dāng)時(shí)，由于可知:由于解的收斂性，則此時(shí)是平庸解。當(dāng)時(shí)，則此時(shí)是平庸解。當(dāng)時(shí)，，其中。所以，，因?yàn)樗?，，則，初始條件：最終，，2數(shù)值解——隱式格式目前，研究熱傳導(dǎo)問題特別是非穩(wěn)態(tài)熱傳導(dǎo)問題十分重要。這里使用隱式格式。利用，關(guān)于t進(jìn)行向前差商：；關(guān)于x進(jìn)行二階中心差：；代入偏微分方程可以得到隱式差分格式：（1）3證明隱式格式的相容性與穩(wěn)定性（1）相容性代入隱式格式得：(2)將（2）與原微分方程相減，得到截?cái)嗾`差：所以此隱式格式與原微分方程相容。（2）穩(wěn)定性令網(wǎng)格比為，則可以將（1）式改寫得到：（3）首先令：（4）將（4）代入（3）式，根據(jù)歐拉公式化簡得：（5）故得放大因子是：所以根據(jù)Fourier方法，隱式格式恒穩(wěn)定。4數(shù)值解——分析與Matlab實(shí)現(xiàn)邊值與初值離散化將邊值與初值離散化，與式（3）聯(lián)立得差分線性方程組：，，再將方程組改寫成的形式：本題的邊界條件均為零。所以可以將上式改寫。Matlab的實(shí)現(xiàn)桿長1米，時(shí)間2秒。設(shè)計(jì)空間步長h=和時(shí)間步長t=，網(wǎng)格比是。從而得到劃分的空間網(wǎng)格點(diǎn)數(shù)是M1+1，時(shí)間網(wǎng)格點(diǎn)數(shù)是M2+1。先設(shè)初始的溫度矩陣U(M2+1,M1+1)。再將邊界條件和初始條件編寫到表示溫度分布的矩陣中。具體代碼可見最后附錄。編寫矩陣A核心代碼：對角線：A(i,i)=1+2r對角線的右方和下方：A(i,i+1)=-r;A(i+1,i)=-r;下面就要運(yùn)用進(jìn)行迭代。當(dāng)k=1時(shí)，A*U(2,j)=U(1,j)當(dāng)k=2時(shí)，A*U(3,j)=U(2,j)當(dāng)k=3時(shí)，A*U(4,j)=U(3,j)以此迭代下去直到k=M2。就可以得到整個(gè)溫度隨時(shí)間和空間的分布矩陣U。數(shù)值解畫圖，如圖1(a)和圖1(b)所示。圖1(a)數(shù)值解的溫度分布圖現(xiàn)在將著色平穩(wěn)過渡。圖1(b)著色平穩(wěn)過渡的數(shù)值解的溫度分布圖5數(shù)值解與解析解的比較首先，我們需要將解析解離散化，解析解中有一項(xiàng)，當(dāng)n越來越大時(shí)，會快速趨于0，故我們可以取n=8000。現(xiàn)在來證明可行性，在matlab里的工作空間運(yùn)算。將解析解的溫度分布畫出來，數(shù)值解畫圖，如圖2所示。圖2解析解的溫度分布圖將數(shù)值解與解析解相減，得到誤差圖。如圖3(a)和圖3(b)，我們從圖3(a)上可以看出空間上的誤差，在邊界處誤差比較大。圖3（a）數(shù)值解與解析解空間誤差我們從圖3(a)上可以看出時(shí)間的誤差，在時(shí)間的最開始，處誤差最大，然后又有一個(gè)小的波動，最后就誤差漸漸變小，最后趨于0。圖3（b）數(shù)值解與解析解時(shí)間誤差6隨時(shí)間變化的細(xì)桿上的溫度分布情況從數(shù)值解的溫度分布三維圖，如圖4(a)和圖4(b)可以看出隨著時(shí)間的增加，細(xì)桿溫度下降最后趨于0℃。從物理角度來說：細(xì)桿的溫度會不斷地向兩端擴(kuò)散，熱量會慢慢散失，最終隨著時(shí)間的增加，細(xì)桿的溫度會趨于0℃。圖4(a)細(xì)桿溫度隨時(shí)間的變化圖現(xiàn)取細(xì)桿中心處一點(diǎn)，觀看它隨時(shí)間的溫度變化情況。圖4(b)細(xì)桿中央（x=）溫度隨時(shí)間的變化圖7穩(wěn)定后細(xì)桿上的溫度分布情況從圖像上可以看出，最后穩(wěn)定的情況下，細(xì)桿的溫度是0℃。參考文獻(xiàn)[1]馮立偉．熱傳導(dǎo)方程幾種差分格式的MATLAB數(shù)值解法的比較[J]．沈陽化工大學(xué)，遼寧沈陽．2011(6)．[2]一維熱傳導(dǎo)方程數(shù)值解法及Matlab實(shí)現(xiàn)[EB/OL]．2014-11-20附錄代碼：%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%此程序用于解決一維熱傳導(dǎo)方程:ut-a^2uxx=0%%邊界條件：u(0,t)=u(L,t)=0%%初始條件：u(x,0)=100,x!=0和L%%u(0,0)=0%%u(L,0)=0%%其中,a^2=1,L=1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%clc;clearall;%區(qū)域及劃分網(wǎng)格L=1;%單位長度的細(xì)桿T=2;%時(shí)間h=;%%%%空間的劃分%%%%t=;%%%%時(shí)間的劃分%%%%r=t/(h*h);%網(wǎng)格比%設(shè)計(jì)步長M1=L/h;M2=T/t;%構(gòu)造邊界條件%構(gòu)造的矩陣：U(時(shí)間，空間）U=zeros(M2+1,M1+1);%編程包含邊值，如U(k,1)=u(0,t)fork=1:M2+1%時(shí)間劃分了M2份，有M2+1個(gè)節(jié)點(diǎn)U(k,1)=0;%兩個(gè)邊界處溫度恒為零U(k,M1+1)=0;end;%構(gòu)造初始條件forj=2:M1%位置劃分了M1份，有M1+1個(gè)節(jié)點(diǎn)U(1,j)=100;end;U(1,1)=0;U(1,M1+1)=0;%差分格式的矩陣形式A*U(k+1,j)=U(k,j)%構(gòu)造矩陣AA=zeros(M1-1);fori=1:M1-1A(i,i)=1+2*r;end;fori=1:M1-2A(i,i+1)=-r;A(i+1,i)=-r;end;%構(gòu)造AU=B中的B%本題邊值的特殊，矩陣B大大簡化了B=zeros(M1-1,1);fork=1:M2j=2:M1;B(j-1,1)=U(k,j);x=A\B;forj=2:M1U(k+1,j)=x(j-1);%k+1時(shí)刻的不同位置的溫度分布end;end;%作圖x=0:h:1;y=0:t:2;[xx,yy]=meshgrid(x,y);figure(1);surf(xx,yy,U);shadingflattitle('一維熱傳導(dǎo)方程--數(shù)值解--溫度分布圖');xlabel('位置x');ylabel('時(shí)間t');zlabel('溫度T');figure(2)s=0;fori=1:8000s=s+(200*(1-(-1)^i))/(i*pi)*sin(i*pi*xx).*exp(-i^2*pi^2*yy);end;surf(xx,yy,s);title('一維熱傳導(dǎo)方程--解析解--溫度分布圖');xlabel('位置x');ylabel('時(shí)間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('一維熱傳導(dǎo)方程--誤差--溫度分布圖');xlab